Collisional formation of detectable exomoons of super-terrestrial exoplanets

Collisional formation of detectable exomoons of super-terrestrial exoplanets

Uri Malamud, Hagai B. Perets, Christoph Schäfer, Christoph Burger
Department of Physics, Technion Israeli Intitute of Technology, Technion City, 3200003 Haifa, Israel
Institut für Astronomie und Astrophysik, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
Department of Astrophysics, University of Vienna, Türkenschanzstraße 17, 1180 Vienna, Austria
Accepted XXX. Received YYY; in original form ZZZ

Exomoons orbiting terrestrial or super-terrestrial exoplanets have not yet been discovered; their possible existence and properties are therefore still an unresolved question. Here we explore the collisional formation of exomoons through giant planetary impacts. We make use of smooth particle hydrodynamical (SPH) collision simulations and survey a large phase-space of terrestrial/super-terrestrial planetary collisions. We characterize the properties of such collisions, including the debris mass, the long-term tidal-stability and the expected properties of the formed exomoons. We find our giant-impact models to generate relatively iron-rich massive disks, however even under the most optimistic assumptions these disks would not form exomoons exceeding Mars mass. Most of the modelled exomoons are likely to form beyond the synchronous radius of the planet, and would tidally evolve outward rather than be destroyed. We also find one rare case in which an exomoon forms through a graze & capture scenario between a super-Earth and an Earth-sized planet. The result is an intact, planet-sized exomoon, containing about half the mass of the Earth. Our results suggest that it is extremely difficult to collisionally form currently-detectable exomoons orbiting super-terrestrial planets, through single giant impacts. It might be possible to form more massive, detectable exomoons through several mergers of smaller exomoons, formed by multiple impacts, however more studies are required in order to reach a conclusion. Given the current observational initiatives, the search should focus on more massive planet categories. However, typical exomoons predicted by our models are very likely to be detectable, given an order of magnitude improvement in the detection capability of future instruments.

planets and satellites: detection, planets and satellites: formation
pubyear: 2019pagerange: Collisional formation of detectable exomoons of super-terrestrial exoplanetsLABEL:lastpage

1 Introduction

For the past two decades, thousands of exoplanets have been identified, providing the first detailed statistical characterization of their properties (Christiansen, 2018). However, to date, there has not been even a single confirmed detection of an exomoon.

The formation of exomoons has been relatively little studied, although they could play an important role in planet formation. Moreover, exomoon environments are important due to their potential of hosting liquid water, thereby creating more opportunities for harbouring life (Williams et al., 1997), and extending the normal boundaries of what is considered habitable environments. Such a possibility is intricately contingent upon multiple factors, including the amount of insolation, tidal heating, and other heat sources available to exomoons (Dobos et al., 2017), as well as their orbital stability (Gong et al., 2013; Hong et al., 2015; Spalding et al., 2016; Alvarado-Montes et al., 2017; Zollinger et al., 2017; Grishin et al., 2018; Hamers et al., 2018; Hong et al., 2018), atmosphere (Lammer et al., 2014; Heller & Barnes, 2015) and the magnetic field of either satellite or planet (Heller & Zuluaga, 2013). Massive exomoons are also important since they can prevent large chaotic variations to their host planet’s obliquity (Sasaki & Barnes, 2014), thereby creating a more stable climate which may be essential to the survival of life on a (solid and in the habitable zone) planet (Nowajewski et al., 2018). It has also been shown that water-bearing exomoons are in principal capable of retaining their water as their host stars go through their high-luminosity stellar evolution phases, and so they can host life or provide the necessary water for supporting life even around evolved compact stars (Malamud & Perets, 2017).

Observationally, several different ways have been proposed for the search of exomoons (Kipping et al., 2009; Simon et al., 2009; Liebig & Wambsganss, 2010; Peters & Turner, 2013; Agol et al., 2015; Noyola et al., 2016; Sengupta & Marley, 2016; Forgan, 2017; Lukić, 2017; Berzosa Molina et al., 2018; Vanderburg et al., 2018). Transit-based techniques are the most promising methods given current observational capabilities, e.g. the Hunt for Exomoons with Kepler (HEK) (Kipping et al., 2012) initiative. The detectability of an exomoon chiefly relies on its orbit, mass and the mass of its host planet (Sartoretti & Schneider, 1999). With the HEK study, exomoons are not likely to be detected below a lower mass limit of about 0.2 Earth masses (). At this mass, any exomoon would be at least one order of magnitude more massive than any satellite in our own Solar system. Such a simple restriction is therefore already suggestive of certain intrinsic properties of the majority of exomoons, given their non-detection so far.

In order to form an exomoon massive enough to be detectable, in accordance with the aforementioned criteria, it is required to have an unusually (by Solar system standards) large satellite-to-planet mass ratio. In addition, one requires a massive host planet, at the very least on the scale of a terrestrial planet, and more likely a Super-Earth or a gas giant. To illustrate the point, an exomoon around an Earth-analogue, requires a satellite to planet mass ratio of 1:5 in order to be detectable, twice than the Pluto-Charon mass ratio, which represents the largest mass ratio known in the Solar system. Typical in-situ formation of satellites inside cirumplanetary disks results in satellite-to-planet mass ratios of the order of (Canup & Ward, 2006), although more recent models (Cilibrasi et al., 2018) show that statistically, more massive satellites can dwell inside the rare tail in the mass distribution. It is therefore an unlikely way of forming detectable exomoons, unless their host planets are in the super-Jupiter mass range. In contrary, giant impacts are readily capable of forming satellites with large satellite-to-planet mass ratios. The most notable examples in the Solar system are the giant impact scenario that formed the Pluto-Charon system (Canup, 2005) and the one that formed the Earth-Moon system (Canup & Asphaug, 2001), with satellite-to-planet mass ratios of and respectively. Such collision geometries involving solid bodies are certainly plausible in the late stages of terrestrial planet formation (Elser et al., 2011; Chambers, 2013), and therefore might credibly give rise to massive exomoons around Earth-like or Super-Earth exoplanets.

Our goal in this paper is therefore to map the collision phase space relevant to the formation of massive exomoons around super terrestrial planets, including new formation pathways which have not yet been suggested in the existing collision formation literature, presented in Section 2. We then briefly introduce in Section 3 the model used for hydrodynamical collision simulations, the considerations for our parameter space, and introduce our pre- and post-processing algorithms. In Section 4 we present the simulation results, and discuss their implications in Section 5, including some predictions of exomoon detections around super-terrestrials, when using present-day or near-future instruments.

2 Collisional formation of exomoons

Most simulation studies of giant impacts have focused on the collisional phase space conductive to the formation of Solar system planets and satellites (Barr, 2016). Despite an extensive collision simulation literature, there have only been a few studies that investigated giant impacts relevant to exoplanets that are more massive than the Earth (Marcus et al., 2010a, b; Liu et al., 2015; Barr & Bruck Syal, 2017). In particular, only Barr & Bruck Syal (2017) (hereafter BB17) focus on the formation of exosolar satellites (or rather the disks from which they accreted), while all the rest examine the effects on the exoplanets themselves.

The study of BB17 examines collisions onto rocky exoplanets up to 10 . Their goal is to identify the collision phase space capable of generating debris disks massive enough to form detectable exomoons by present-day technology. They use a Eulerian Adaptive Mesh Refinement (AMR) CTH shock physics code, and examine different masses, impact geometries and velocities. While BB17 manage to demonstrate, for the first time, detailed simulations capable of forming very massive proto-satellite disks – only 2 cases in a suite of 28 impact simulations (7%) result in a disk mass of 0.3 , hence beyond the 0.2 criteria. Furthermore, the disk mass is only a hard upper limit on the mass of the exomoon that could form. The actual fraction of mass in the disk that coagulates into a moon is not immediately clear. It primarily depends on the initial specific angular momentum of the disk, and also on how one chooses to model the disk.

In N-body disk models that assume a particulate disk of condensed, solid particles, thus neglecting the presence of vapour, about 10–55% of the mass of the disk would go into the satellite (Kokubo et al., 2000). In more complex hybrid models that consist of a fluid model for the disk inside the Roche limit and an N-body code to describe accretion outside the Roche limit, about 20-–45% of the mass of the disk would go into the satellite (Salmon & Canup, 2012). The hybrid models provide a more realistic view since gravitational tidal forces both prevent accretion inside the Roche limit and also disrupt planet-bound eccentric clumps, scattered from outside the Roche limit by close encounters. A sufficiently energetic giant impact dictates that this zone is initially in a state of a two-phase liquid/vapour silicate ’foam’ and as such its evolution is controlled by the balance between viscous heat dissipation (further inducing vaporization) and radiative cooling (Ward, 2012). This disk spreads inward toward the planet, and outward beyond the edge of the Roche limit, where newly spawned clumps and their corresponding angular momentum join the particulate, satellite accreting disk. The complexity of such models, however, entails large uncertainties on the disk physics (Charnoz & Michaut, 2015). The problem is more accentuated when one considers the conditions applicable to the formation of detectable exomoons around super-terrestrial planets, in which more massive or hotter disks are involved. Ignoring such complications, that is, if the aforementioned accretion studies are scalable to larger masses, it would suggest no more than about half the mass of the disk would form the final satellite. Given that assumption, the two outlier cases in the study of BB17 form exomoons in which the final mass is less than 0.15 , hence below the detecion criteria.

The primary goal of this study is therefore not only to reproduce, but also extend the BB17 collision phase space in an attempt to identify more likely host planets or paths to forming massive exomoons. In the following paragraphs we consider: (a) more giant impact scenarios which we speculate are likely to form massive disks; and (b) the formation of massive exomoons through consecutive, multiple impacts, rather than in a single giant impact.

In order to include additional collision simulations that might increase the chance of forming massive exomoons, we look more closely at the BB17 collision phase space for specific hints. Figure 1 is a reproduction of Figure 6 in the study of BB17, and shows their results for the disk-to-total mass ratio () as a function of the normalized angular momentum (). The latter is given by Canup (2005):


where is the impact angle, the impact velocity, the escape velocity, the impactor-to-total mass ratio and . As can be seen, and are well-correlated, with the exception of six outlier cases that result in very low disk masses and do not fit the rest of the data. We have identified the outlier cases in the study of BB17 to be consecutively, velocity6 to angle3 in Table 1. The latter correspond to all the cases with either a high velocity or a low impact angle. In other words – high velocity and low impact angle appear counter conductive to the formation of massive debris disk, and thus we judge them as incompatible avenues to forming massive exomoons.

Figure 1: A reproduction of Figure 6 from BB17, plotting as a function of . 6 outlier cases in a suite of 28 collisions produce low mass disks and do not fit with the rest of the data, having either high impact velocities or low impact angles.

In the upper part of Table 1, 19 simulations have approximately the canonical Lunar-forming value (0.11). The remaining 9 (in a total of 28 simulations) have much larger values, and they form 7 of the 9 disks with the highest disk-to-total mass ratios. We thus judge large to generate favourable outcomes, perhaps unsurprisingly, as they are more compatible with the Pluto-Charon impact scenario. In Section 3.2 we describe our additional simulations with a large value of .

As previously mentioned, we also consider the stepwise growth of exomoons through multiple, rather than a single impact. This sequence of impacts is simply part of the critical collisional evolution that naturally takes place during the last stages of terrestrial-type planet accretion. Multiple impacts create multiple exomoons. The tidal evolution and migration of these formed exomoons and their mutual gravitational perturbations (Rufu et al., 2017; Citron et al., 2018) could then result in several possible evolutionary outcomes, with roughly equal probabilities, including collisions among two exomoons (eventually growing into a more massive final exomoon), ejection of exomoons, or their recollisions with the host planet (Malamud et al., 2018). Estimating the mass of such an exomoon, formed through mergers, is beyond the scope of this paper, since it requires a large set of N-body and moon-moon impact simulations. Assuming that such a moon can form, however, subsequent impacts onto the planet could in principal clear the moon by either ejecting it, or triggering a moonfall. For the latter case, we wish to investigate what would be the likely result. Based on the statistical analysis of Citron et al. (2018), most moonfalls have extremely grazing geometries, with 90°, as shown in Figure 2. For an Earth-sized planet, Malamud et al. (2018) have shown that extremely grazing moonfalls do not always result in the complete-loss of moon material, but instead they often give rise to a new generation of intact moons that retain most of their original mass. If this behaviour is shown here to be typical of super-terrestrial planets as well, then in the framework of multiple-impact formation, moonfalls may have a high probability of continuing the process of exomoon growth, rather than restarting it.

Figure 2: The distrbution of moonfall impact angles triggered by the gravitational perturbations between an inner and an outer moon (=0°: head-on impact; =90°: exteremely grazing impact) based on analysis of data from the Citron et al. (2018) study.

Finally, we note that in the course of this study, we will have an opportunity to perform a comparison between impact simulations using our Lagrangian SPH method, and the Eulerian AMR method in the BB17 study. This in itself has some important value, since (to our knowledge) only one other comparison study has been performed until now for planetary impacts (Canup et al., 2013), and it involved only Lunar-forming scenarios. While a fully consistent comparison is difficult when the precise details involved (set-up procedures, pre-processing and post-processing algorithms) are handled by different groups (see e.g., Section 3.6), we can nevertheless show if using the SPH or AMR methods broadly result in the same overall trends in the data, and in particular the same amount of exomoons or disks which are potentially detectable using present-day instruments.

3 Methods

3.1 Hydrodynamical code outline

We perform hydrodynamical collision simulations using an SPH code developed by Schäfer et al. (2016). The code is implemented via CUDA, and runs on graphics processing units (GPU), with a substantial improvement in computation time, on the order of several times faster for a single GPU compared to a single CPU, depending on the precise GPU architecture. The code has already been successfully applied to several studies involving impacts (Dvorak et al., 2015; Maindl et al., 2015; Haghighipour et al., 2016; Wandel et al., 2017; Schäfer et al., 2017; Burger et al., 2018; Haghighipour et al., 2018; Malamud et al., 2018).

The code implements a Barnes-Hut tree that allows for treatment of self-gravity, as well as gas, fluid, elastic, and plastic solid bodies, including a failure model for brittle materials. Given the analysis of Burger & Schäfer (2017) however, and the typical mass and velocity considered for our impactors and targets (see Section 3.2), we perform our simulations in full hydrodynamic mode, i.e., neglecting solid-body physics, being less computationally expensive. We use the M-ANEOS equation of state (EOS), in compatibility with the BB17 study. Our M-ANEOS parameter input files are derived from Melosh (2007).

3.2 Parameter space

We are exploring the phase space of potential large-exomoon forming impacts, which are more likely to be detectable, even with more advanced instruments than currently available.

As a starting point, we are repeating the impacts considered in the study of BB17, albeit using a different numerical method. All other parameters being equal, we can determine whether using an SPH code instead of the AMR code of BB17, might in itself improve or impede the formation of massive exomoons. We therefore analyse the BB17 suite of collisions, using identical composition and EOS. The BB17 suite of collisions is listed in the upper part of Table 1.

As discussed in Section 2, we also extend the BB17 collision phase space. We identified giant-impact scenarios with a high impactor-to-total mass ratio, as being conductive to forming massive disks. Accordingly we set 4 new simulations similar to pc_earth and big_earth in the study of BB17, albeit with more massive planets, and also with equal mass impactor-target scenarios, as seen in the lower section of Table 1 (bigger1-2 and biggest1-2).

To complete our suite of simulations, we consider six additional moonfall cases, whereby existing exomoons with masses 0.05, 0.1, 0.2, 0.3, 0.4 and 0.5 (moonfall1-6) respectively, are gravitationally perturbed to collide with their host super-Earth planet at 90°, as discussed in Section 2 (in the context of multiple-impact exomoon growth). For simplicity we assume that their impact velocity equals the mutual escape velocity , being a reasonable upper limit, according to the analysis of the impact velocity distribution data from the (Citron et al., 2018) study.

3.3 Initial setup

We consider differentiated impactors and targets composed of 30% iron and 70% dunite by mass. Both impactors and targets are non-rotating prior to impact, and are initially placed at touching distance. Including rotation as a free parameter would entail a very large increase in the number of simulations, especially if one assumes a non-zero angle between the collisional and equatorial planes. Here we assume no rotation and a coplanar collision geometry, in order to comply with the BB17 study. As we discuss in Section 4.2, the initial rotation in particular could have some significance, and we suggest to explore it in future dedicated studies (see Section 5.3).

The initial setup of each simulation is calculated via a pre-processing step, in which both impactor and target are generated with relaxed internal structures, i.e. having hydrostatic density profiles and internal energy values from adiabatic compression, following the algorithm provided in appendix A of Burger et al. (2018). This self-consistent semi-analytical calculation (i.e., using the same constituent physical relations as in the SPH model) equivalently replaces the otherwise necessary and far slower process of simulating each body in isolation for several hours, letting its particles settle into a hydrostatically equilibrated state prior to the collision (as done e.g., in the work of Canup et al. (2013) or Schäfer et al. (2016)).

Altogether we have 38 simulations, using a resolution of SPH particles. The simulations were performed on the bwForCluster BinAC, at Tübingen University. The GPU model used is NVIDIA Tesla K80. Each simulation ran on a single dedicated GPU for approximately 7-10 days on average, tracking the first 28 hours post collision to comply with BB17 on what they considered a complete simulation. In total we thus have 54 weeks of GPU time.

3.4 Debris-disk mass analysis

Planet scale impact simulations often result in a cloud of particles, some of which will accrete onto the planet, other remain in a bound proto-satellite disk, and some escape the system. In order to determine which particles go where, an analysis is performed as a post-processing step. Our algorithm resembles the one used by BB17. BB17 refer to the procedures described by Canup et al. (2001) as the basis for their analysis. Given the latter, we note that some of the details required in order to compare our approaches are missing or are incompatible with our interpretation of the data (explanation will be provided below). As a consequence, we point out that there may be minor variations in our respective analyses results, however we judge them to have a small effect because the overall scheme follows essentially a similar classification approach.

Our detailed classification algorithm follows this 5-step algorithm:

(a) We find physical fragments (clumps) of spatially connected SPH particles using a friends-of-friends algorithm.

(b) The fragments are then sorted in descending order according to their mass.

(c) We classify these fragments in two categories: gravitationally bound (GB) to the planet and gravitationally unbound (GUB). The first fragment (i.e., the most massive, in this case the target/proto-planet) is initially the only one marked as GB, and the rest are marked as GUB.

(d) We calculate:


where and are the center of mass position and velocity of GB fragments, denoting indices of GB fragments and , and are the corresponding mass, position and velocity of each fragment.

Then, for each fragment marked as GUB we check if the kinetic energy is lower than the potential energy:


where is the gravitational constant and is the summed mass of gravitationally bound fragments. , and are the fragment mass, position and velocity, denoting indices of GUB fragments. If equation 3 is satisfied then fragment is switched from GUB to GB. We iterate on step (d) until convergence (no change in fragment classification within an iteration) is achieved. At this point we deem the mass of the planet as equal to that of most massive GB fragment, while the other GB fragments will be classified as either belonging to the planet or to the bound disk according to the following, final step.

(e) We calculate the respective fragment energy , eccentricity and pericentre , denoting the indices of GB fragments, however excluding P (planet bound) fragments, i.e., deemed as belonging to the planet:


and are the P fragments’ centre of mass coordinates and are calculated in a similar way to equation 2. Since we assume that the collision geometry is coplanar to the planet’s equator, each fragment that satisfies , i.e., has its pericentre residing inside the planet’s physical cross-section – may be classified as P, having a planet-bound trajectory. We iterate on step (e) until convergence is achieved. The calculation of is specified in Section 3.5.

Once the algorithm is complete we have the unbound mass , planet-bound mass and disk mass , respectively. We also track for each one, their respective compositions and the origin of SPH particles (impactor or target).

3.5 Tidal stability analysis

In addition to calculating the mass of the disk (to check its potential of forming a massive exomoon), BB17 also discuss if a satellite can survive its post-accretion tidal evolution. The Roche radius is the distance below which tidal forces would break the satellite apart, and therefore the initial accretion distance of the satellite has to be greater. The Roche radius is given by , where is the planet’s radius, the planet’s density and the satellite’s density. The mean distance of satellite formation in giant-impact simulations is slightly beyond the Roche radius, at around 1.3 (Elser et al., 2011).

After the initial accretion, the satellite’s orbit evolves by tidal interactions with the planet, depending on its position relative to the synchronous radius , the distance at which the satellite’s circular orbital period equals the rotation period of the planet, or in other words, the satellite’s orbital mean motion () equals the initial rotation rate of the planet (). The synchronous radius depends on the rotation rate of the planet, and can be obtained by equating the gravitational acceleration () and the centripetal acceleration ( or ), giving .

Satellites that form inside the synchronous radius, orbit their planets faster than their planet rotates. In this case the tidal bulge they raise on the planet lags behind, and acts to decelerate them, and so they spiral towards the planet and are ultimately destroyed when they cross the Roche radius. However, satellites that form outside the synchronous radius, recede from the planet as angular momentum is tidally transferred from the planet to the satellite. In either case, more massive satellites tidally evolve faster. One can reach a conclusion on the tidal fate of a satellite, comparing the two radii. Since in giant-impacts a satellite typically forms from a disk just outside the Roche radius, it implies that if , the satellite evolves outwards and survives. It is also understood that in slowly rotating planets, moves outwards, therefore making it much more difficult to form tidally-stable satellites. We note that the aforementioned discussion applies to prograde satellites. Retrograde satellites which orbit in the opposite sense relative to the planet’s rotation, will of course also in-spiral. In this paper however, we start all our simulations with non-rotating targets and impactors, hence all satellites will orbit in the same sense.

To complete the set of equations, the initial rotation rate of the planet can be calculated from the angular momentum of the material judged to constitute its mass. For this calculation we take only the mass of the largest fragment . We note that by the end of the simulation, the small fraction of planet-bound material which has not yet accreted, cannot contribute more than a few % to the final planet mass anyway, and is not expected to change significantly. For we then find the total angular momentum by the summation of its individual particle angular momenta, where denotes SPH particle mass and and denote SPH particle relative (to fragment’s center of mass) position and velocity. We then calculate the angular momentum unit vector , which points towards the direction of the rotation axis. In order to get the rotation rate we calculate , the distance vector from the rotation axis to the particle relative position. Then the total moment of inertia is similarly given by , and the rotation rate by .

For a rotating object, we expect to have some degree of flattening , such that:


where and are its equatorial and polar radii respectively. For such an object with mass , the density is generally given by:


In order to calculate the planet’s density we thus need to know its equatorial and polar radii. However the latter two quantities cannot be calculated directly before all planet-bound material has accreted, so we will assume for simplicity that . Such an assumption is judicious since the mass of the planet is, as previously mentioned, very close to the mass of the biggest fragment at the end of our simulations, while non-negligible changes to the density typically entail a large change in the mass (hence the pressure by self-gravity). By the same consideration, and since changes in were also regarded as negligible, we will assume consistently that the flattening . Now can be directly calculated from Equation 8 substituting for . The equatorial and polar radii of the biggest fragment are physically computed by considering its constituent particles, such that and sml (the smoothing length distance). From the latter we can calculate using Equation 7. We note that these equatorial and polar radii values, which are measured relative to the planet rotation axis, are simply indicative of the planet’s shape, and therefore not always identical to those of a hydrostatically flattened ellipsoid (see Figure 3 and associated discussion). In that sense is simply the ratio given by Equation 7, which helps us calculate the planet’s physical cross-section.

In contrast to the planet, the satellite has not fully accreted at the end of the simulation, because the time scale for full accretion is at least two orders of magnitude longer. We thus have no direct information about its physical properties, but we do have information about the disk from which it will coagulate. In order to calculate we will use the fact that the satellite is not nearly as massive as the planet, and therefore as a heuristic approach, it could be more readily estimated from , denoting the indices of the disk’s constituent materials, the relative material mass fractions and their corresponding specific densities. We note that unlike in Solar system satellites, for extremely massive exomoons we expect the actual to be somewhat larger than this estimation, and therefore our ensuing should be considered an upper limit.

Finally, we rewrite Equations 7 and 8 for and :


The effective planet radius is then:


Equations 9-11 are iteratively re-evaluated in step (e) of Section 3.4, in order to obtain the disk mass.

3.6 Algorithm differences from previous studies

Since we will be comparing our results with previous studies, we wish to note as a caveat that the technique for calculating the equatorial radius in the BB17 study, and therefore the planet and disk masses, is based on an algorithm from an Earth-related study by Canup et al. (2001), wherein an iterative process is used to estimate these radii from the Earth’s mass and density, assuming the latter is equal to the known present-day Earth density value. This value is of course unsuitable for the much more massive exo-planets considered by BB17, and would lead to a ’too-small’ planet radius, and therefore Roche radius, which also bears directly on the tidal stability conclusion. BB17 however did not specify what density value they did use in their study, nor the details of how it was calculated. Our calculation of the planet’s density is however based on measuring it directly from its physical size and mass, as was described above.

Another difference in the algorithm is in the calculation of the planet’s flattening. Canup et al. (2001) assume, as we do, that during the initial disk formation process, only a clump of material whose pericentre resides inside the equatorial radius, will merge with the planet (see step (e) of Section 3.4). However, their approach is to calculate this physical cross-section from the planet’s flattening coefficient, which is in turn calculated according to a prescribed formula (Kaula, 1968). The latter is given as a function of the planet’s rotation rate, which can be calculated from its angular momentum. Our inspection of the data, however, shows that for extremely energetic impacts the actual planet’s post-collision shape is not always that of a standard flattened ellipsoid in hydrostatic equilibrium (unlike in canonical Lunar-forming giant impacts). In some cases, we see an extended region of relatively high density material (>1000 kg m) at temperatures of several 10000 K, as in Figure 3. Previous studies argue that the relevant viscous timescale of such material is much longer than the initial disk formation timescale, dominated by gravity (Ward, 2012), hence the effective cross-section for clump mergers can be larger relative to the one obtained by using the prescribed, rotation-dependent flattening coefficient formula from Canup et al. (2001). We therefore define the flattening coefficient according to Equation 7, where the planet’s equatorial and polar radii are measured with respect to its data-derived shape, as previously described.

Given the arguments above, we caution that there may be some differences primarily between our planet radius calculation and that of BB17. Therefore the Roche radius and to a lesser extent also the disk mass analysis, may result in somewhat different values, despite using mostly identical equations and procedures.

Figure 3: Edge-on view of the planet at the end of the gamma3 simulation. The (semi-transparent) colour scheme shows the planet’s density in kg m. The inner planet is surrounded by a cloud of relatively high-density (>1000 kg m), and extremely hot silicate material (up to 40000 K). The planet’s collision cross-section does not comply with a standard flattened ellipsoid shape.

We also point out that a few of our collisions do not necessarily form a disk of debris, but rather result in graze & capture scenarios, which form intact moons. In that case, the Roche radius is unrelated to tidal stability because it is not necessarily indicative of the initial distance of the satellite. It is rather the satellite’s pericentre that we must compare to , in order to gain some approximative understanding of how the orbit will develop. Such considerations will be discussed in the following Section.

4 Results

The results from all simulations are summarized in Table 1, which shows, from left to right, the parameters for different scenarios (upper table, reproduction of the BB17 suite of collisions and lower table, new collisions in this study); the disk mass in the BB17 study; and the disk mass + synchronous & Roche radii + composition in this study.

Parameters BB17 This study
big_earth 2.311 0.325 0.929 62.13 0.127 0.068 1.162 2.087 0.586 0.357
pc_earth 1.013 0.327 0.983 61.92 0.061 0.084 0.962 1.478 0.755 0.465
rocky_exo2 0.266 0.146 0.986 50.55 0.006 0.003 0.805 1.07 0.428 0.195
rocky_exo3 0.466 0.142 0.962 50.55 0.006 0.007 0.938 1.29 0.429 0.209
rocky_exo4 2.324 0.123 0.862 50.55 0.04 0.033 1.669 2.203 0.435 0.302
rocky_exo5 3.805 0.119 0.825 50.55 0.032 0.053 1.915 2.587 0.45 0.279
rocky_exo6 7.208 0.115 0.774 50.55 0.033 0.03 2.165 3.301 0.34 0.298
rocky_exo7 11.151 0.111 0.736 50.55 0.058 0.011 2.439 3.973 0.169 0.458
rocky_exo8 18.066 0.106 0.691 50.55 0.063 0 2.851 4.83 0 0.522
ser119 1.015 0.130 0.922 50.55 0.014 0.008 1.246 1.7 0.369 0.328
velocity3 7.208 0.115 0.864 50.55 0.112 0.098 2.307 3.233 0.41 0.276
velocity4 7.208 0.115 0.953 50.55 0.071 0.046 2.253 3.289 0.337 0.355
velocity5 7.208 0.115 1.042 50.55 0.311 0.506 3.515 2.937 0.678 0.103
velocity6 7.208 0.115 1.216 50.55 0.003 0.001 4.34 3.073 0.498 0.274
velocity7 7.208 0.115 1.390 50.55 0.005 0 5.27 3.256 0.249 0.402
velocity8 7.208 0.115 1.647 50.55 0.002 0 6.615 3.181 0.337 0.344
angle1 7.208 0.115 0.774 20.70 0.001 0.001 3.438 3.555 0 0.448
angle2 7.208 0.115 0.774 28.11 0 0.002 2.875 3.554 0 0.276
angle3 7.208 0.115 0.774 36.09 0.001 0.004 2.503 3.553 0 0.22
angle4 7.208 0.115 0.774 62.07 0.09 0.154 2.322 3.178 0.469 0.197
angle5 7.208 0.115 0.774 70.46 0.151 0.163 2.278 3.162 0.486 0.281
gamma1 7.126 0.196 0.778 50.55 0.112 0.092 1.775 3.302 0.31 0.296
gamma2 7.131 0.321 0.778 50.55 0.12 0.087 1.571 3.396 0.179 0.529
gamma3 7.264 0.451 0.771 50.55 0.298 0.14 1.563 3.402 0.197 0.559
gamma4 7.518 0.442 0.757 50.55 0.136 0.137 1.572 3.418 0.23 0.541
gamma5 7.176 0.137 0.775 50.55 0.063 0.054 2.025 3.293 0.339 0.306
gamma6 7.112 0.258 0.779 50.55 0.135 0.106 1.645 3.411 0.158 0.33
gamma7 7.182 0.386 0.775 50.55 0.148 0.135 1.561 3.382 0.209 0.471
bigger1 7.2 0.327 0.983 61.92 0.286 1.565 3.009 0.622 0.08
bigger2 7.2 0.5 0.983 61.92 0.196 1.719 3.485 0.054 0.462
biggest1 11 0.327 0.983 61.92 0.262 1.715 3.833 0.253 0.22
biggest2 11 0.5 0.983 61.92 0.299 1.967 4.03 0.015 0.511
moonfall1 7.2 0.007 1 89 0.012 120.6 2.986 0.718 0
moonfall2 7.2 0.014 1 89 0.014 50.5 2.994 0.685 0
moonfall3 7.2 0.028 1 89 0.19 43.41 2.941 0.724 0
moonfall4 7.2 0.042 1 89 0.241 15.92 2.913 0.745 0
moonfall5 7.2 0.056 1 89 0.331 14.94 2.899 0.747 0
moonfall6 7.2 0.069 1 89 0.431 14.52 2.888 0.744 0

The upper part of the table consists of the 28 collisions performed by BB17 and repeated in this study. The lower part of the table lists the extended suite of collision simulations considered only in this study. Columns from left to right: (1) simulation name ; (2) total mass of target and impactor in Earth mass units ; (3) impactor-to-target mass ratio ; (4) impact velocity in units of the mutual escape velocity ; (5) impact angle ; results from BB17: (6) disk mass in Earth mass units ; results from this study: (7) disk mass in Earth mass units ; (8) synchronous rotation radius ; (9) Roche radius ; (10) disk iron fraction ; (11) disk target-material fraction. For more details on the notations see Section 3.2. The highlighted disk mass values mark the cases in which it nears or exceeds the 0.2 detection criteria.

Table 1: Suite of collisions.

4.1 Disk and exomoon masses

Results from this study show that only one disk from the original BB17 suite of 28 collisions, velocity5, has a mass larger than the 0.2M detection criteria. A close inspection of the data however, shows this to be a unique case whose outcome is a graze & capture scenario (See Figure 4). All other simulations in the upper part of Table 1 generate a disk of debris, and the final exomoon mass is subject to further considerations/modelling (see Section 2 for discussion). For velocity5 the disk mass equals the mass of the intact exomoon. Any additional accretion from the disk onto the exomoon is negligible, however at 0.506 M (and 61% of the mass of the original impactor), the captured exomoon is already well above the HEK detection criteria.

Most of our new scenarios in the lower part of Table 1, generate disk masses that are approximately 0.2 M or more. As we shall discuss in the following Section (4.2), the bigger1-2 and biggest1-2 disks give rise to tidally stable exomoons, however the detectability of such moons with HEK is still in question, since even in the best case scenario, at least 2/3 of the mass of these disks must coagulate to form a moon in order for it to exceed the HEK detection criteria. For terrestrial-scale collisions, no more than about 1/2 of the mass of the disk would coagulate to form a moon (see Section 2). Further studies are thus required for super-terrestrial-scale collisions (see Section 5.3).

The moonfall simulations have mixed results. Analysis of the data shows moonfall1-2 to be graze & run scenarios. These moons graze the planet, however there is insufficient dissipation to capture them, and the most massive fragments are unbound. In contrast, moonfall3-6 are graze & capture scenarios, and retain most of the mass of the original impacting moons (95%, 80%, 83% and 86%, respectively).

(a) 0.17 hours
(b) 1.5 hours
(c) 3 hours
Figure 4: First 3 hours of the velocity5 collision scenario: an original 0.83 M impactor, forms a detectable 0.5 exomoon through a graze & capture scenario. The (semi-transparent) colour scheme shows density in SI units. Resolution is particles.

4.2 Tidal stability

BB17 do not provide the detailed tidal stability analysis results for all their simulations. They mention only that pc_earth is marginally stable against planetary tides with and that velocity5 has . In our results pc_earth is also stable to planetary tides, with , whereas velocity5 has . The difference between the studies could arise from the calculation of , as mentioned in Section 3.6. Moreover, since the velocity5 exomoon forms in a graze & capture scenario, rather than coagulates from a disk, the Roche radius does not imply where is the initial location of the satellite. Therefore stability has to be evaluated based on its orbit, which can be directly extracted from the simulation data. We did not follow the entire post-impact orbit of this exomoon with SPH for computational-cost reasons, however its pericentre can be calculated based on its trajectory and it was found to be 3.392 km. The latter is larger than its , so it would not disrupt tidally. It is also of the order of , which implies long-term stability.

We find that all other impacts in the upper part of Table 1, apart from velocity6-8, are stable to tidal interactions. Whereas velocity6-8 are hit-and-run scenarios, and the angular momentum carried away by the debris increases with the velocity, resulting in a slower planet rotation and therefore moves outwards, inhibiting stability. Had the planet been rotating, already prior to the impact, perhaps high-velocity collisions would have been able to produce tidally stable moons.

We also note that velocity5 and moonfall3-6, while having a different collision angle, are all graze & capture scenarios. They thus similarly create intact exomoons with relatively little debris, such that nearly all the mass is concentrated in the largest fragment. We are able to analyse the data in order to calculate their pericentre distances : 33924 km, 51874 km, 69230 km, 137922 km and 146211 km, respectively. The latter are all found to be larger than their respective . In other words, none of these exomoons will experience a breakup by tidal forces on their second pass. The pericentre distances of velocity5 as well as moonfall5-6 are also of the order of , or larger, which implies long-term stability. In contrast, moonfall3-4 do not fulfil the stability requirement because their is too large for their pericentre distances. However, we must keep in mind that the large is a result of the slow post-impact rotation of the initially non-rotating planet, generated by the relatively small-sized impactors. In reality, the planet is far more likely to be rotating prior to the moonfall, which will bring further inwards. For example, pre-rotation of the moonfall3 and moonfall4 targets with periods of 12h and 19h respectively, would be sufficient to make their captured exomoons tidally stable, even without considering the extra angular momentum transferred by the moonfall. For comparison, note that nearly all the planets in the upper part of Table 1 have lower than 51874 km, and thus have a shorter rotation period than 12h, which strongly suggests this is not an uncommon state.

4.3 Composition

Our results indicate that disks generated by the energetic impacts considered here for super-terrestrial planets, are often (a) iron-rich and (b) are composed mostly from impactor materials, although in about a third of the cases the impactor/target material fractions are almost identical. The former result is rather dissimilar to impacts around Earth-like planets, which noticeably generate extremely iron-poor disks, independent of whether SPH or AMR methods are used (Canup et al., 2013). Impactor material fractions are however not so unlike.

Information about the origin of disk material is not provided in the study of BB17. The iron disk fraction is only mention once, for velocity5. Without providing an exact fraction, they broadly suggest that the velocity5 exomoon is iron-rich, which validates our result. In fact, in our entire suite of simulations the captured moon from velocity5 is among the richest in iron. With nearly 70% iron, it has a similar composition to planet Mercury, which also may have had much of its mantle stripped by a single impact or multiple impacts (Benz et al., 1988; Chau et al., 2018), however it is approximately nine times more massive than Mercury. The other graze & capture scenarios, moonfall3-6, are similarly iron-rich.

5 Discussion

5.1 Comparison with previous studies

This work follows only one previous study which was recently carried out by BB17, and examined the possibility of forming exomoons detectable by HEK. Together, these are the only two studies ever performed which consider the formation of exomoons, given the highly energetic impacts relevant to super-terrestrial planets. In Section 3.6 we argue that the details concerning the post-processing analysis, may lead to differences in the results between the two studies, which will not be related to the numerical method used, namely, SPH versus AMR respectively. We nevertheless observe that the disk masses follow precisely the same trend and that they vary by up to a factor of two, and typically much less, as can be directly compared in Table 1. The differences are more noticeable, and also more important for the most massive disks, because it is only the latter that have the potential to form exomoons which are potentially detectable using current technology.

The most prominent example involves scenarios gamma3 and gamma4. The latter have very similar collision parameters, and are expected to yield similar outcomes, as they do, in fact, in our study. In the study of BB17, however, the gamma3 disk is twice as massive as that of gamma4. In Figure 3 which shows a snapshot of the density profile for the gamma3 post-impact planet, we discuss how the disk mass calculation is sensitive to how one models the planet’s cross-section. Normally, we use a different analysis calculation compared to BB17 since we derive the appropriate scale-length directly from the data (due to the unorthodox shapes that sometimes emerge), rather than from a prescribed flattening (see Section 3.6). If we nevertheless change our algorithm to take the planet radius based on its material-derived density (in the same way we calculate it for the exomoons), our analysis outcome becomes a lot more similar to that of BB17. Hence, we see that the simulation method itself is sometimes not as important as the details of analysis, and even then differences between the two studies never amounted to more than a factor of two. A more accurate comparison requires (a) intimate knowledge of analysis details; (b) comparison between multiple parameters – not only disk mass; and (c) visual inspection of data. It is therefore hard to manage if not performed by the same authors, in addition to the fact that only our study offers full details on the composition, origin and stability of each exomoon/disk. Overall, we nevertheless conclude that the two methods broadly result in similar outcomes, somewhat validating the more accurate SPH versus AMR comparison study made by Canup et al. (2013).

5.2 Detectable exomoons: predictions

The main goal of this paper was to predict the likelihood of forming, and therefore finding observable exomoons, now or in the future. Given the simplest HEK detection restriction, it requires an exomoon which is at least 0.2 . A previous study by BB17 set a more modest goal of identifying the collision phase space capable of generating debris disks more massive than 0.1 . They found that only 10 cases in a suite of 28 simulations were able to generate disks over 0.1 , and only 2 cases generated disks over 0.2 , the latter being the actual restrictive condition from Kipping et al. (2009). BB17 acknowledge that the disk mass can only serve as an upper limit on the final exomoon mass, and yet they conclude that "detectable rocky exomoons can be produced for a variety of impact conditions". As shown in Table 1, we have been able to both reproduce the BB17 results with similar outcomes, and also greatly increase the number of disks with masses over 0.2 . Yet based on their results, as well as on ours, we actually arrive at exactly the opposite conclusion, that is, with respect to the HEK criteria.

In Section 2 we show how studies of moon accretion restrict the fraction of mass in the disk that coagulates into a moon, to about 10-55%. These studies were not performed for massive disks around super-Earths, and yet the chance to accrete an exomoon detectable by HEK from a disk which is only slightly more massive than 0.2 appears marginal at best, even for the most optimistic assumptions. In this study we have hypothesized, based on the results of BB17, that the most likely avenues to generate massive disks are through Pluto-Charon type giant-impacts, scaled up to the super-Earth size range. We have tested this hypothesis and it turned out to indeed exclusively generate massive disks over 0.2 . It is however not entirely clear how plausible or frequent such collisions are. They probably represent only a tiny fraction of the potential collision phase space during terrestrial planet formation. Moreover, since these disks never amount to more than 0.3 , the final mass of the exomoons generated by these disks should be in the range 0.03-0.15 . Hence, below the HEK detection criteria.

We therefore conclude that only one giant-impact simulation, velocity5, resulted in an exomoon whose mass is well above the HEK detection criteria. This scenario involves a collision between a 7.2 Super-Earth and a 0.83 , Earth-like planet. We have analysed the data and found that the exomoon forms intact, through a graze & capture collision, keeping most of the original impactor mass (0.5 ). We find this moon to be marginally stable to tides and conclude that it will most likely survive its long-term dynamical evolution around its host-planet. Overall, the relatively large characteristic impactor-to-total mass ratio, the typical collision angle of 50° and the collision velocity at approximately , all require very specific conditions. While different combinations and small variations in mass, collision angle and velocity, will likely result in a similar outcome, we nevertheless conclude that the velocity5 captured exomoon must represent the exception, rather than the rule.

We also consider, for the first time, the possible formation of an exomoon detectable by HEK through several mergers between smaller exomoons, instead of in a single giant-impact. Multiple impacts are a natural consequence of late terrestrial planet accretion, and have been considered previously as a possible formation mechanism for the Earth’s moon (Citron et al., 2018). Each impact forms a smaller moon which tidally evolves and gravitationally interacts with an existing, previously formed moon. With roughly equal probabilities, these moons can eject, merge, or fall back onto the planet. For the latter case, they typically fall at extremely grazing geometries with impact angle of the order of 90° (Malamud et al., 2018). Here we check whether such moonfalls lead to the destruction, ejection or survival of the impacting exomoons. We find that moonfalls of impactors below approximately 0.2 become unbound, while massive impactors in excess of 0.2 are likely to be captured and retain most of their mass. In other words, an exomoon already above this mass limit can survive a moonfall, increasing its probability to continue its growth in the subsequent planet impact. However, the only way to form it to begin with through several mergers of smaller (and prevalent) exomoons, requires that the gravitational interaction between them would exclusively result in mergers, while avoiding ejections and moonfalls, since both would terminate its growth. Further studies are required in order to establish the relative probabilities of mergers, ejections and moonfalls, since the latter depend on the relative masses of the interacting moons, and since studies so far (Citron et al., 2018) have targeted the Earth rather than super-Earth planets. With targeted studies and statistics, we may be able to reach a decisive conclusion.

According to this reasoning, the best and only way of forming a sufficiently massive exomoon to be detectable by HEK, requires a rare graze & capture scenario involving a super-Earth and an Earth-like planet. In the coming decade, only PLATO might enable a slight improvement, being able to detect Mars-like, 0.1 exomoons (private correspondence with David Kipping). Unless we considerably improve our ability to detect exomoons with future instruments, our results indicate predict a great difficulty in finding the first exomoon around any terrestrial or super-terrestrial planet. On the other hand, if we improve the detection threshold by one order of magnitude to approximately 0.01 , we argue that exomoons are already extremely likely to be detected, since about a third of the disks formed in Table 1 are more massive than 0.1 . The exomoons which will coagulate from these disks will thus be in the mass range 0.01-0.05 (Section 2), hence, detectable to such instruments. Our results indicate these exomoons might be more iron-rich compared to solar system moons.

Our study therefore motivates the development of a new generation of instruments, and sets a specific goal for its developers. Meanwhile, we suggest focusing our efforts with Kepler data on the biggest known Super-Earths, or else more massive planet categories. For the latter, future studies of exomoon formation are definitely required.

5.3 Future studies

Our main suggestions for future studies are as follows:

  • Pre-rotation and collision geometries - both the BB17 study and our study consider initially non-rotating targets or impactors, as well as coplanar impact geometries, for simplicity. It is clear that initial rotation and an inclined collision plane would affect the disk masses and angular momenta. With improved high performance computing capabilities, future studies may choose to account for these complexities and consider a larger grid for the potential collision phase space.

  • Disk evolution and exomoon mass estimation - both the BB17 study and our study result in a wide range of disk conditions, including some very massive disks with extreme temperatures and pressures. It is unclear how such disks evolve to form exomoons and what would be their final masses. As discussed in Section 2, there are several approaches to this problem, some of which are very complex. A relatively simple way to obtain lower-limit estimations of exomoon masses is to self-consistently hand over SPH particles beyond the Roche to be simulated for long-term dynamical evolution in N-body codes. We are currently developing the appropriate tools necessary, and plan to perform such studies in the future.

  • Dynamical and tidal evolution of multiple exomoons - as mentioned in Section 5.2 (in the context of exomoon growth through multiple impacts) we require a detailed statistical understanding of gravitational interactions among exomoons, including their tidal evolution about their Super-terrestrial host planets. Such a study can easily follow the design of Citron et al. (2018), with the required modifications.

  • Exomoon formation around gas-giant planets - our results suggest that exomoons may be more readily detected around more massive categories of planets. There is already one promising candidate, a Neptune-sized exomoon orbiting a super-Jupiter planet (Teachey & Kipping, 2018). The work of Hamers et al. (2018) analytically calculates a tidal capture and subsequent orbital evolution scenario, however until now there are no studies performed to identify the detailed collision-formation of exomoons around this class of planets.

6 Acknowledgment

UM and HBP acknowledge support from the Minerva center for life under extreme planetary conditions, the Israeli Science and Technology ministry Ilan Ramon grant and the ISF I-CORE grant 1829/12. CS and CB acknowledge support by the high performance and cloud computing group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHP and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG. CB acknowledges support by the FWF Austrian Science Fund project S11603-N16.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description