# Survival of Exomoons Around Exoplanets

V. Dobos, S. Charnoz, A. Pál, A. Roque-Bernard, Gy. M. Szabó

► To cite this version:

V. Dobos, S. Charnoz, A. Pál, A. Roque-Bernard, Gy. M. Szabó. Survival of Exomoons Around Exoplanets. Publications of the Astronomical Society of the Pacific, 2021, 133, 10 pp. 10.1088/1538-3873/abfe04. insu-03590015

HAL Id: insu-03590015

https://insu.hal.science/insu-03590015

Submitted on 26 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Survival of Exomoons Around Exoplanets

V. Dobos¹,²,³, S. Charnoz⁴, A. Pál²,⁵, A. Roque-Bernard⁴, and Gy. M. Szabó³,⁵

¹Kapteyn Astronomical Institute, University of Groningen, 9747 AD, Landleven 12, Groningen, The Netherlands; vera.dobos@rug.nl

²Konkoly Thege Miklós Astronomical Institute, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), 1121, Konkoly
Thege Miklós út 15-17, Budapest, Hungary

³MTA-ELTE Exoplanet Research Group, 9700, Szent Imre h. u. 112, Szombathely, Hungary

⁴Université de Paris, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France

⁵ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, Szombathely, Szent Imre h. u. 112, Hungary

Received 2021 February 3; accepted 2021 May 4; published 2021 August 31

## Abstract

Despite numerous attempts, no exomoon has firmly been confirmed to date. New missions like CHEOPS aim to characterize previously detected exoplanets and potentially discover exomoons. In order to optimize search strategies, we need to determine those planets which are the most likely to host moons. We investigate the tidal evolution of hypothetical moon orbits in systems consisting of a star, one planet, and one test moon. We study a few specific cases with ten billion years integration time where the evolution of moon orbits follows one of these three scenarios: (1) “locking,” in which the moon has a stable orbit on a long timescale (≥10⁹ yr); (2) “escape scenario” where the moon leaves the planet’s gravitational domain; and (3) “disruption scenario,” in which the moon migrates inwards until it reaches the Roche lobe and becomes disrupted by strong tidal forces. Applying the model to real cases from an exoplanet catalog, we study the long-term stability of moon orbits around known exoplanets. We calculate the survival rate which is the fraction of the investigated cases when the moon survived around the planet for the full integration time (which is the age of the star, or if not known, then the age of the Sun). The most important factor determining the long-term survival of an exomoon is the orbital period of the planet. For the majority of the close-in planets (<10 days orbital periods) there is no stable orbit for moons. Between 10 and 300 days we find a transition in survival rate from about zero to 70%. Our results give a possible explanation for the lack of successful exomoon discoveries for close-in planets. Tidal instability causes moons to escape or being tidally disrupted around close-in planets which are mostly favored by current detection techniques.

Unified Astronomy Thesaurus concepts: Exoplanets (498); Exoplanet dynamics (490); Tidal interaction (1699)

## 1. Introduction

Regular moons are predicted to form generally in planetary systems, as a direct outcome of planet formation. During the core-accretion, gravitational perturbations between planet embryos imply a series of constructive impacts up to the formation of a fully grown planet (Nagasawa et al. 2007). This phase may happen between a few Myr to a few 100 Myr after the star formation (see e.g., Chambers & Wetherill 1998), and in these processes, satellites may form around the growing planets. Such giant impacts are advocated for the formation of the Earth’s Moon (Canup 2004) and for the formation of Uranus and Neptune’s satellites (Morbidelli et al. 2012). Satellites may also be formed in rings, which could be natural outcomes of giant impacts or tidal disruptions, either in the planet formation phase or later in a

relaxed planetary system (see e.g., Charnoz et al. 2009; Canup 2012; Crida & Charnoz 2012).

An alternative scenario invokes the accretion of moons in the gaseous circumplanetary envelopes that surround the most massive giant planets during their growth in the gaseous protoplanetary disks (see e.g., Canup & Ward 2006; Sasaki et al. 2010; Miguel & Ida 2016). Inside the circumplanetary disk, the solid material is replenished by the surrounding protoplanetary disk. This solid material is thought to accrete in the form of large satellites in a way similar to planets. Then the young satellites migrate inward and sometimes sink into the planet. When the circumplanetary disk disappears, the satellite system remains and then undergoes long-term tidal evolution.

After the planets have been formed, the orbits of moons evolve due to the planet’s tides. Tidal forces lead to the dissipation of energy. On the one hand, this usually results in the gradual expansion of the moon’s orbits, if the evolution starts outside the planet’s synchronous radius. (A synchronous orbit is where the moon’s orbital period is equal to the planet’s rotational

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
period.) On the other hand, moons inside the synchronous orbit will migrate inwards and can reach the planet's Roche limit where the moons are disrupted by tidal forces.

The timescale of orbital evolution depends mainly on the intensity of dissipation in the planet's interior, which relies on the material composition and internal structure of the planet. The dissipation coefficient also varies with the changing planet radius (for gas-rich planets Recent papers attempting to compute the quality factor, $Q$, (also known as *tidal dissipation*) show that $Q$ could even be frequency dependent, and thus may strongly depend on the orbital distance (see e.g., Ogilvie & Lin 2004). However, in the current stage of our knowledge, the large uncertainties in the origin of the dissipation process and in the internal structure of detected exoplanets preclude any theoretical calculation of $Q$ for a specific planet. For this reason, in the present paper, a classical approach is adopted where $Q$ will be considered as a constant, and many values will be explored. The dissipation scales as $Q^{-1}$ and we can assume it to be in a steady-state when exploring those scenarios which lead to moons on the long-term stability. For rocky planets, $Q$ is low and the dissipation is high ($Q \sim 10-100$), conversely for giant planets $Q$ is high ($\sim 10^4-10^5$) and dissipation is low (Goldreich & Soter 1966; Alvarado-Montes et al. 2017).

The quest for exomoon discovery is a challenging task in itself, and upon finding an exomoon we could empirically constrain the $Q$ value of the host planet by considering the planetary distance and the age of the star. This will provide invaluable information about the exoplanet interior, too, since high/low dissipation rates may depend on the planet core and interface with mantle (Lainey et al. 2012; Remus et al. 2012).

Recently, Alvarado-Montes et al. (2017) investigated the orbital changes of exomoons around close-in giant planets due to tidal effects. In addition to the tidal evolution, they implemented the changes of the giant planet's physical parameters into the model. The radius of the planet decreases over time, and this change is most prominent in the first billion years after the formation of the planet. The value of $k_2/Q$ slowly increases, which also has an impact on the orbital evolution of the moon.

Probably the most intriguing question concerning the moons is the lack of a firmly confirmed moon in exoplanetary systems. Before novel observation projects are designed for exomoon detection, first we need to find the answer to the question of where the exomoons would be. Did we use suboptimal signal detection strategies, or did we simply look for moons in systems where their existence is not too likely? This latter possibility is the central question of this paper.

Barnes & O'Brien (2002) showed that exomoons are more stable at longer planetary orbital periods, and also that the orbit of exomoons can tidally evolve first outward and then inward

(toward the planet) as the planetary rotation slows due to tidal interaction with the star. Heller (2012) studied the orbital stability of habitable moons using the method of Kipping (2009) and it was found that for stars of masses below $\approx 0.2 M_\odot$, moons are not likely to stay on stable orbits around planets which are in the habitable zone. Later Zollinger et al. (2017) showed that for planets in the habitable zones of $\lesssim 0.5 M_\odot$ stars, the orbit of moons is not stable as they migrate toward the planet.

Here we investigate the evolution of hypothetical satellites orbiting the currently known exoplanets with numerical integration, under the influence of tidal interaction. We map a wide range of parameters, including planet periods, moon to planet mass ratios, and $k_2/Q$. Instead of studying a few cases in large detail and with high accuracy, which would be meaningless due to our lack of knowledge on both the physical origin of tidal dissipation and on the planet's internal structure, we are more interested in the statistical survival rate of moons. In order to investigate a large range of configurations over billions of years of evolution, we used the average description of Barnes & O'Brien (2002) which are formulated in Equations (1)-(3) (Section 2). We apply this method to identify different evolutionary pathways of moons (Section 3) and then we calculate survival rates of hypothetical moons around known exoplanets with the aim of pointing out specific planets as observation targets that might have satellites (Section 4 and the supplementary data).

Independently of our work, Sucerquia et al. (2019) used similar methods to investigate the orbit stability of exomoons around close-in planets, and estimated their detectability through TDV and TTV effects. Guimaraes & Valio (2018) had a similar goal with a focus on the photometric detectability of hypothetical moons around Kepler planets. They selected 54 planets that could possibly host observable moons. Although we also list planets with survival rates for moons, we focus our investigations on survival statistics.

## 2. Methods

### 2.1. Orbital Evolution of Satellites

We calculate the orbital evolution of the planet and the moon from the following equations (Barnes & O'Brien 2002; Sasaki et al. 2012).

$$ \frac{dn_m}{dt} = -\frac{9 k_{2p} R_p^5 G M_m n_m^{16/3}}{2 Q_p (GM_p)^{8/3}} \operatorname{sgn}(\Omega_p - n_m), \quad (1) $$

$$ \frac{dn_p}{dt} = -\frac{9 k_{2p} R_p^5}{2 Q_p G M_p (GM_*)^{2/3}} n_p^{16/3} \operatorname{sgn}(\Omega_p - n_p), \quad (2) $$
$$ \frac{d\Omega_p}{dt} = -\frac{3}{2} \frac{k_{2p} R_p^3}{Q_p} \frac{(GM_m)^2}{\alpha(GM_p)^3} n_m^4 \operatorname{sgn}(\Omega_p - n_m) \\ -\frac{3}{2} \frac{k_{2p} R_p^3}{Q_p} \frac{1}{\alpha GM_p} n_p^4 \operatorname{sgn}(\Omega_p - n_p), \quad (3) $$

where $n_m$ and $n_p$ are the mean motions of the moon and the planet, $k_{2p}$ is the second-order Love number of the planet which relates to the dilatation due to tidal response, $Q_p$ is the quality factor of the planet which describes the dissipation of energy in the body, $G$ is the gravitational constant, $R_p$ is the radius of the planet, $M_m$, $M_p$ and $M_*$ are the masses of the moon, the planet, and the host star, respectively, $\Omega_p$ is the rotation frequency of the planet, and $\alpha$ is the prefactor of the moment of inertia ($\alpha M_p R^2$). Circular orbits are assumed for both the planets and the satellites.

Since Equations (1) and (2) have the form of $\dot{y} = C y^p$ (where $y$ is either $n_m$ or $n_p$, $p$ is the exponent of $16/3$ while $C$ contains all of the astrophysical constants) and Equation (3) can directly be integrated once the time evolution of $n_m$ and $n_p$ are known, this set of equations can be solved analytically on the domains until $n_m$, $n_p$ or $\Omega_p$ do not surpass each other; in this case, the presence of the sgn() function would not allow the fully analytical solution. We, therefore, implemented the solution in a piecewise-analytical manner, monitoring the sign changes in the differences between $n_m$, $n_p$ or $\Omega_p$ while in the distinct pieces, the exact solution is evaluated. Thus, besides the vicinity of the points where the signs commute, the solution does not even depend on this stepsize, allowing us to highly speed up the computations.

For calculating the survival rate of moons around known exoplanets, we let the orbits of the planet and the moon evolve until the current age of the star, or if it is not known then until the age of the Sun (4.57 Gyr). We solve the equations analytically with fixed time steps to monitor whether the moon reaches the Roche lobe of the planet or the critical distance (half of the Hill sphere). In these cases, the integration ends because the moon is considered to be lost either because of tidal disruption or escape from the planet's gravitational bond (for the latter case see Domingos et al. 2006). The critical distance is calculated by $R_{crit} = 0.5 R_{Hill} = 0.5 a_p (M_p/(3M_*))^{1/3}$.

If at the end of the calculation the moon is still orbiting the planet, then it is counted as a “survival” case.

## 2.2. Parameters used for the Integrations

Stellar mass ($M_*$), planet mass ($M_p$), semimajor axis ($a_p$) and radius of the planet ($R_p$) were taken from The Extrasolar Planets Encyclopaedia (http://exoplanet.eu, data obtained on 2021 March 25). Only those planets are considered where the $M_p$ or the $M_p \cdot \sin i$ or the $R_p$ value is known (at least one of them). Planets with higher masses than 13 Jupiter mass are ignored as well as those with $M_* < 0.08 M_\odot$. In addition, the selected planets also had

to have either the orbital period ($P_p$) or the semimajor axis data. After this selection 4064 planets remained.

For both the radius and mass (or $M_p \cdot \sin i$ if $M_p$ is not given) of the planet, the uncertainties were taken into account wherever possible. If there was no uncertainty in the data, then simply the given mass or radius was used. If there was a symmetric uncertainty given, then the mass or radius values were generated by choosing a random value with a Gauss distribution. The symmetry was determined with the same method as described in Chen & Kipping (2017, Section 2.2) where the difference between the upper and lower errors must be below 10%. For these symmetric cases, the mean of these two values was taken as the standard deviation of the data. If the difference was higher than 10%, then the errors were ignored and simply the mean value was used. We applied exactly the same method for those planets where only the $M_p \cdot \sin i$ value was known, not the real mass. In these cases, the error of $M_p \cdot \sin i$ was taken into account wherever possible. This way we could generate realistic data for all planets, taking into account the uncertainties in the measurements.

There are 1007 cases (out of 4064) where both the radius and the mass ($M_p$ or $M_p \cdot \sin i$) of the planet are known. For all the other planets, the Forecaster model was used for generating realistic values for the missing parameter (Chen & Kipping 2017). This tool predicts the mass (or radius) of the planet based on the radius (or mass). As an input, the generated (above-mentioned Gauss distributed) mass (or radius) values were used as posteriors. This way, for each run the planet had a different mass or radius, or both (if the original data had symmetric errors).

See the following description of the defined ranges and distribution of the chosen random values for the rest of the (initial) parameters.

1. **Mass of the moon:** The mass of the hypothetical moons is selected randomly with a uniform distribution between 1% and 10% of the host planet's mass, for each planet.

2. **Semimajor axis of the moon:** The initial semimajor axis of the moon is chosen randomly with a uniform distribution on a logarithmic scale between $2 R_p$ and $R_{crit}$.

3. **Radius of the moon:** Estimating the radius of exomoons is very challenging because only the solar system serves as an example. We choose to define it through their densities ($\rho_m$) using their mass which is already chosen. Based on the big moons (radius >200 km) in the solar system, we define three categories depending on the locations with respect to the snowline ($a_{snow}$). This is because beyond the snowline there are icier, hence less dense bodies in the solar system. The snowline is calculated as $a_{snow} = T_{eff}^2 \cdot R_* / T_0^2$ where $T_{eff}$ is the effective temperature of the star and $T_0 \approx 230$ K is the equilibrium temperature at the planets sub-stellar point (Cowan & Agol 2011). A fourth category is also given for those
cases when $T_{\text{eff}}$ or $R_*$ is not known, hence $a_{\text{snow}}$ cannot be calculated. The $\rho_m$ values are selected randomly with a Gauss distribution with the following parameters.

(a) If $a_p < a_{\text{snow}}$: the mean value of $\rho_m$ is $3 \text{ g cm}^{-3}$, $\sigma = 1/3 \text{ g cm}^{-3}$.

(b) If $a_{\text{snow}} \le a_p < 2a_{\text{snow}}$: the mean value of $\rho_m$ is $2.5 \text{ g cm}^{-3}$, $\sigma = 1/3 \text{ g cm}^{-3}$.

(c) If $2a_{\text{snow}} \le a_p$: the mean value of $\rho_m$ is $2.5 \text{ g cm}^{-3}$, $\sigma = 1/6 \text{ g cm}^{-3}$.

(d) If $a_{\text{snow}}$ is not known: the mean value of $\rho_m$ is $2.5 \text{ g cm}^{-3}$, $\sigma = 1/2 \text{ g cm}^{-3}$.

4. *Rotation period of the planet:* In the solar system the planet that has the slowest spin around its axis is Jupiter with 9.9 hr. Scholz et al. (2018) define an empirical spin-mass relation for planets and brown dwarfs. According to their work, the spin velocity scales as $v \sim \sqrt{M}$, and from this relation it can be concluded that the spin period does not exceed 5 days even for brown dwarfs. (Spin-orbit resonances are not taken into account.) Based on this information and considering possible deviations, we set the spin period to a random value with a uniform distribution between 10 hr and 5 days.

5. *Quality factor of the planet:* For the quality factor ($Q_p$) three categories are defined.

(a) *Rocky planets* if $R_p < 2R_⊕$: $10 < Q_p < 500$ (Goldreich & Soter 1966). Since lower values are more probable, the distribution of the randomly chosen $Q_p$ value is uniform on a logarithmic scale.

(b) *Ice/gas giants* if $R_p \ge 2R_⊕$ and $P_p > 10$ days: $10^3 < Q_p < 10^6$ (Goldreich & Soter 1966; Ogilvie & Lin 2004; Lainey et al. 2012, 2017). For these planets intermediate values are the most probable in this range, hence $Q_p$ values are chosen with a log-normal distribution with a mean of $Q_p = 10^{4.5}$ and $\sigma = 10^{0.5}$.

(c) *Hot Jupiters* if $R_p \ge 2R_⊕$ and $P_p \le 10$ days: $Q_p \sim 5 \times 10^6$ (Ogilvie & Lin 2004). We do not have much information on the possible $Q_p$ value of these planets, hence we use a Gauss distribution with $Q_p = 5 \times 10^6$ as the mean value and $\sigma = 2 \times 10^6$.

$Q_p$ is a critical parameter, as it controls the tidal evolution of the planet, however, its value is very poorly constrained, and, in the current state of our knowledge, it could vary by several orders of magnitudes even for similar planets (Remus et al. 2012). For the giant planets in the solar system, $Q$ is mostly constrained by the secular acceleration of a big moon orbiting relatively close to the planet and also by theoretical considerations. For these planets, $Q$ is estimated to be in the range of $10^4-10^5$ (see e.g., Goldreich & Soter 1966; Ogilvie & Lin 2004; Lainey et al. 2009). However, if $Q$ is not treated as a constant but as a function of the tidal frequency, then we obtain much lower values for Saturn with $Q$ about 2500 (Lainey et al. 2012, 2017). $Q$ is mostly unknown for ice giants like Uranus

and Neptune, but theoretical considerations put them in the range $>10^4$ (Goldreich & Soter 1966). However, this estimate is based on some unconstrained assumptions on the moon formation scenarios and ages. There is also a more recent estimation for Neptune which puts the $Q$ value between 9000 and 36,000 (Zhang & Hamilton 2008).

Conversely, the $Q$ of terrestrial planets is very low, in comparison to gas giants, and is typically between ten and a few hundreds (see e.g., MacDonald 1964; Goldreich & Soter 1966; Smith & Born 1976; Ray et al. 2001; Rainey & Aharonson 2006; Lainey et al. 2007). The $Q$ for these planets is better constrained because of the possibility of laboratory tests of rocks (where possible). While $Q$ might go up to 500 for rocky planets, it is typically below 200 in the solar system.

It is highly challenging to draw conclusions about the possible $Q$ value of exoplanets. It seems reasonable to assume that rocky planets have similar values as in the solar system. We can also speculate that maybe the $Q$ of giant planets is in the same order of magnitude as their solar system counterparts, however, theoretically, with low-density ice-rich cores, $Q$ in the range of a few hundreds would also be possible for a Saturn-like planet (Lainey et al. 2017). In hot Jupiters, because of the assumed synchronous spin and circularized orbit, $Q$ might be around $5 \times 10^6$ (Ogilvie & Lin 2004).

The Love number $k_2$ varies modestly from one planet to another in the solar system but is between 0.1 and 0.6 (Gavrilov & Zharkov 1977). We choose $k_{2p} = 0.5$. However, $k_2$ never appears alone, but only in $k_2/Q$. So we fix the value of $k_2$ and vary the value of $Q$ to save a free degree of liberty that would lengthen the computation time.

For the same reason, the parameter used for calculating the gyration radius is also fixed to $\alpha = 0.3$. It is expected that planets with different compositions and structures have different gyration radii, but by keeping $\alpha$ fixed, a free parameter is removed.

For each planet, ten thousand runs are made to get a wide sample of the different initial parameters. What we call the “survival rate” is simply the ratio of surviving moons (those that stayed in orbit around the planet until the end of the integration) and the total number of tested configurations.

For a few example cases, the parameters of the planet and the moon are set to specific values which are described in detail in Section 3.

### 3. Example Cases

In the example cases (see Figures 1 and 2), the system configuration consists of a Solar-mass star with an Earth-like planet and a moon that has a different mass in the different cases. The quality factor of the planet is $Q_p = 10$ in all cases. The initial parameters used for the different runs are summarized in Table 1. Note that the time shown in the horizontal axis is significantly shorter in Figure 1 (where it ends at $5 \times 10^6$ yr) than in Figure 2 (where it ends at $10^{10}$ yr).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Moon orbit and planet spin evolution in two example cases. See the initial parameters in Table 1. The Blue dashed horizontal line at $6 \times 10^{-6} s^{-1}$ represents the location of the critical distance from the planet (half of the Hill radius of the planet) beyond which the moon escapes. The Red dashed horizontal line shows the Roche radius of the planet expressed in the unit of mean motion. Upper panel (disruption scenario): the moon quickly synchronizes with the planet, but when migrating inwards to the planet, it reaches the Roche lobe and becomes tidally disrupted. Lower panel (escape scenario): the moon slowly migrates outwards meanwhile the planet spin evolves. When the moon reaches the critical distance (after about 4 million years), it escapes from the planet.</figcaption></figure>

### 3.1. Evolution Outcomes

In the integrations, the moon has been observed to follow three different families of evolutionary tracks, and depending on the actual parameter set, all of these three tracks can evolve rapidly, resulting in an unstable system, or the evolution can be slow enough to lead to a surviving moon. Both the morphology and the stability depend on the position of the Hill sphere, thus eventually, on the masses and orbital periods of both the planet and the moon, the stellar mass, and $Q_p$.

#### 3.1.1. Disruption Scenario

When the planet is close enough to the star, the tidal torque on the planet from the star very effectively spins down the planet. The tidal bulge is shifted slightly toward the trailing side of the moon's orbit, the moon quickly loses energy, spirals inward, and is disrupted at the Roche lobe. The timescale of this migration is very short (in the order of 1–200 million years,

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Moon orbit and planet spin evolution in three example cases. See the initial parameters in Table 1. Blue and red dashed horizontal lines represent the location of the critical distance for moon escape and the location of the Roche limit, respectively, expressed in the unit of mean motion. Upper panel (moon–planet locking scenario): the moon synchronises with the planet after $\sim 3 \times 10^9$ yr; up to this point, the moon slowly migrates outwards from the planet, but after synchronisation the moon’s orbital period slowly decreases. Middle panel (star–planet locking scenario): for a smaller $M_m/M_p$ mass ratio the planet synchronises with the star and becomes tidally locked. Bottom panel (temporarily locked scenario): first the planet synchronises with the star (at around $2 \times 10^9$ yr) but because of the moon’s inward migration, later it synchronises with the moon instead ($\sim 7 \times 10^9$ yr), but after that, the moon shortly reaches the Roche limit of the planet.</figcaption></figure>

mostly depending on the stellar distance and on the moon–planet mass ratio) compared to the 10 Gyr integration time, therefore we consider that these moons are unstable. This
<table><caption>**Table 1**
Initial Parameters Used for the Example Cases Shown in the Different Panels of Figures 1 and 2</caption><thead><tr><th></th><th>Figure 1<br>Upper<br>Panel</th><th>Figure 1<br>Lower<br>Panel</th><th>Figure 2<br>Upper<br>Panel</th><th>Figure 2<br>Middle<br>Panel</th><th>Figure 2<br>Bottom<br>Panel</th></tr></thead><tbody><tr><td><i>M</i><sub>m</sub>/<i>M</i><sub>p</sub></td><td>1/85</td><td>1/85</td><td>1/85</td><td>1/150</td><td>1/300</td></tr><tr><td><i>a</i><sub>p</sub> (au)</td><td>0.3</td><td>0.3</td><td>1.0</td><td>1.0</td><td>0.8</td></tr><tr><td><i>n</i><sub>m</sub> (10<sup>-6</sup><br>s<sup>-1</sup>)</td><td>38.0</td><td>6.1851</td><td>2.6640</td><td>2.6640</td><td>2.6640</td></tr></tbody></table>

Note. The star mass is $M_{*} = 1M_{\odot}$, the planet mass is $M_{p} = 1M_{\oplus}$ and the quality factor of the planet is $Q_{p} = 10$ in each case.

scenario is described in detail in Alvarado-Montes et al. (2017), and our results are completely compatible with their findings (Figure 1, top panel).

### 3.1.2. Escape Scenario

If the planet's orbit is close to the star, then the Hill radius is smaller. On any orbit inside the Hill sphere, the moon suffers strong tidal forces, especially if the planet is relatively massive, which makes the outward migration very fast. As a consequence, if the moon formed close to the critical distance from the planet (half of the Hill sphere), then it will quickly escape before synchronization could happen (Figure 1, lower panel).

### 3.1.3. Tidal Locking

In connection with tidal locking, we define three different scenarios.

**Moon-planet locking: the Selenoid scenario.** The evolution track is somewhat similar to the Earth-Moon case. First, the moon migrates outwards, until a spin-orbit synchronisation is reached between the orbital period of the moon and the rotation of the planet. In the later evolution, the moon slowly turns back, the typical timescale before colliding with the planet was in the order of 10 billion years in our integrations. (More specifically, the end state of the moon is reaching the Roche radius of the planet, where the moon will be disrupted). Therefore, we consider this scenario as stable (Figure 2, top panel). We note that the similarity with the Earth-Moon case is limited to the synchronization between the two bodies, but there is a qualitative difference: here the planet spin is synchronizing with the orbit of the moon, while in the Earth-Moon case the spin of the Moon is synchronized with the orbit of the Earth.

**Star-planet locking.** This scenario is similar to the Selenoid, but rather than synchronising with the moon, the planet is synchronised to the central star. We did find this possibility in several test runs and seems to be typical for certain configurations (close-in systems with the lightweight moon). In these cases, the moon's migration is slow and has a similar timescale as in the Selenoid case. The path is qualitatively identical: first,

the moon migrates outwards and then later inward (Figure 2, middle panel). This configuration is also stable in the long term, and the moon can be observed (with difficulties coming from its low mass and small size).

**Temporarily locked scenario.** We found another tidal evolution scenario with a rather complex evolution for which we have not found an earlier discussion in the literature. In the beginning, it has essentially the same evolution as the star-planet synchronisation, but since the moon migrates inwards, here the tidal forces from the moon is increasing, and from a point, they majorate over the tidal forces from the star. At this point, the planet rotation will be decoupled from the mean motion and resynchronizes with the moon. The inward migration of the moon is enhanced from that point, and the moon soon reaches the Roche limit of the planet. In the scenario shown at the bottom panel of Figure 2, the planet's synchronisation with the star happens after approximately $2 \times 10^9$ yr, where there is the turning point of the moon's orbit, too, and the break-up of the star-planet synchronisation occurs at $\sim 7 \times 10^9$ yr. Because of the large timescales, this configuration is also stable in the long term, and moons in such systems could be found.

### 3.2. Timescales

The timescales of the discussed scenarios depend on the migration rate, hence indirectly on the distance of the Hill-sphere which can limit the evolution path very close to the star.

Close-in planets always evolve very rapidly, because (i) the Hill-sphere is close to the planet, and the moon reaches the critical distance from the planet very soon or turns back at a very close distance to the planet; and (ii) since the moon is always very close to the planet, very prominent tidal forces emerge and the orbit evolution will be very fast. For most close-in configurations, similarly to the hottest sub-Jupiters known, the tidal escape timescale of the moon can be even less than 1000 yr. The only possible scenarios for those cases when the planet is too close to the star are the disruption and the escape scenarios (Figure 1 upper and lower panels).

If the planet orbits the star at a larger distance, the path to the critical distance from the planet is longer and the evolution is slower, hence the timescale of stability can reach billions of years, or even the Hubble timescale (these are the tidal locking scenarios, see Figure 2). For selenoids, this means that the fallback takes several billion years, and for runaway moons, the runaway time takes several billion years. In other words, with increasing planet orbital period first we get quickly falling-back selenoids, then stable selenoids (Figures 1 and 2 upper panels).

## 4. Known Exoplanets

We calculated the evolution of possible model moons in the known planetary systems with the aim of pre-selecting those
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The survival rate of moons around known exoplanets as a function of the planet's orbital period. Black colored dots: planets with known radius and mass data. Teal colored dots: planets without known radius data, for these, their radius was estimated based on their measured mass (or $M_p \sin i$). Yellow-colored dots: planets without known mass data, for these the masses were estimated from their measured radius.</figcaption></figure>

exoplanets which can physically have a satellite. These planets might worth a detailed follow-up observation for a moon survey. Planets and their model moons were simulated according to the recipe described in Section 2. The results are listed in Table A.1 in the supplementary data, available online at stacks.iop.org/PASP/133/094401/mmedia, ordered by decreasing survival rate. (In the table only those cases which have non-zero survival rates are shown.) The table also contains the mean of the integration times for each planet. If the survival rate is a hundred percent, then the mean of the integration times will be equal to the age of the star (or of the Sun if the star's age is not known, as described in Section 2). If in some cases the integration stops earlier, then the mean of the integration times will be shorter. This parameter together with the survival rate can be considered as a measure for the stability of moon orbits.

Figure 3 shows the survival rate for test exomoons as a function of the orbital period of the planet. Because the critical escaping distance is small below 10 days orbital period, only a few exomoons stay in a stable orbit around the planet. Above ~300 days, the survival rate is about 70%-90%, as a consequence of the large critical distance. Between about 5 and 300 days we find a "transition zone" where the increasing size of the critical escaping distance results in higher and higher survival rates.

The saturation at ~85% survival rate for large orbital periods is caused by the chosen initial semimajor axes of the moons which in some cases is below the synchronous orbit. In these cases, the moons spiral inwards to the planet until they reach the Roche sphere.

Only a fraction of known exoplanets are promising exomoon hosts according to this work. This is primarily because the orbital period behaves as a moon "censor": not many satellites are expected (survival rate mostly below 40%) for planets in shorter period orbits than 10 days. This is not surprising, since planets, this close to the star have small Hill radii, hence stable orbits for a moon are limited to a very narrow space. As the orbit of the moon evolves, sooner or later it will reach the instability limit and escape from the planet's gravitational bond.

Interestingly, the planets with known mass and radius data (see the black dots in Figure 3) are showing smaller deviation above ~100 days orbital period, compared to those for which the radii or masses are estimated from their other parameter (teal and yellow colored dots, respectively). For the planets with known mass and radius value, the survival rate reaches 50% at around 100 day orbital period of the planet, while for the other planets it is much more diffused in between approximately a hundred and a thousand days. Choosing a mass or radius value for the planet is an extra free parameter that naturally causes higher uncertainty in the results. This might be the reason for the shift of the *transition zone* from 5–100 days to 30–300 days orbital period for the two kinds of planets. For those where we have to add random mass or radius values, we will find more unstable cases, hence the survival rate starts increasing only at higher orbital periods where the critical escape distance from the planet is larger. Another explanation for the shift could be an observation bias, since the planets for which both mass and radius data are measured are those for which both transit and radial velocity observations were made. This means that they are orbiting close to the star and in the
meantime, they are also massive. This might also explain why the transition zone is closer to the star for these planets.

For known exoplanets with long orbital periods, the host stars are quite faint, and the S/N ratio required for a successful moon detection is not achievable by even as precise instruments like the Kepler telescope. This finding is the resolution of the seeming paradox of missing exomoons in Kepler data. Moons are more likely to have a stable orbit around long orbital period planets, but these cases are difficult to observe. This result is also an alert for later transit finding surveys to focus on planets around bright stars with $\gtrsim 10–100$ days orbital period.

## 5. Discussion

The successful detection of a moon around an exoplanet relies basically on two independent factors: the moons must exist (occurrence), and those moons must be observable. The occurrence is a product of two factors: the formation rate and the stability timescale. Formation studies of moons (and planets or stars) are beyond the scope of this paper, however, with the described formalism, we estimate the survival rate of moons around currently known planets.

The survival rate of model moons in known exoplanetary systems shows that satellites can only stay in orbit around the host planet if the orbital period of the planet is larger than five days. In our investigations, we found a transition zone between 5 and 300 days orbital periods where the survival rate is growing with larger orbital periods. Above 300 days orbital period the moon survival rate reaches 70%–90% for the investigated planets. These specific numbers however depend on the range of chosen initial parameters. Also, increasing the integration time could result in decreasing the survival rates in some cases. This would be only a quantitative change as it would not alter the relation between the survival rate of moons and the orbital period of the planets. Our investigation serves statistical purposes, and it shows a clear trend in survival rates that increase with the planetary orbital period.

This result is in good agreement with the findings of Alvarado-Montes et al. (2017), who developed a more detailed model including the changes of the planet’s physical parameters over time as an additional factor influencing the orbital evolution of satellites. Their results showed that moons of planets with lower than ~60 day orbital periods collide with the planet, but for planets with larger stellar distances, the satellites slowly migrate outwards and stay inside the Hill-sphere of the planet during the whole 4 Gyr integration time (“realistic” case, Figure 9 in their article). Including the change in the planet radius and $k_2/Q$ tidal parameter over time in their model increased the chance of moon survival. This implies that neglecting these factors in our study might underestimate the moon survival rates.

Our results are in line with those of Sucerquia et al. (2019), too, who concluded that long-period planets are more favorable

moon hosts than short-period ones. They used a very similar calculation method with the aim of estimating the TTV and TDV signals for detection. They conclude that in their investigated phase space no system produces TDV signals strong enough to be observed, but there are very few cases where the TTV signal would be observable with TESS and Kepler. Larger moon-to-planet ratios and stellar distances could help observations.

With a similar aim, but using a different method, Guimara'es & Valio (2018) also investigated known planets which could host moons. Besides checking the possible semimajor-axes of moons, they also took into account their detectability, providing a list of 54 Kepler planets and candidates which could have detectable moons.

Recently, Tokadjian & Piro (2020) made a similar work to ours in which they investigated the orbit evolution of moons around known exoplanets using simplified tidal lag models. They identified 36 habitable zone planets around which a moon might be stable for long timescales (at least for 1 Gyr). We found that from these planets 29 have at least a 50% survival rate for exomoons (typically between 50 and 65%). There are also four planets with survival rates higher than 65%: these are Kepler-459 b (69%), Kepler-456 b (72%), Kepler-1654 b (70%) and Kepler-1647 b (69%). They also found that for 21 planets the timescale before losing the moon is larger than the Hubble time in all investigated cases. Our results show moon survival rates typically between 55% and 70% for these planets with mean integration times ($t_{mean}$) between 2 and 3 Gyr (see the supplementary data). There is one outstanding case (Kepler-1628 b) that resulted in only a 9% survival rate in our study. The explanation of our lower stability levels is that a much bigger phase space was explored for the initial parameters compared to the study by Tokadjian & Piro (2020). Also, since the age of these planets is typically lower than 5 Gyrs (and in some cases, it is not known), our calculation stopped much earlier.

Martínez-Rodríguez et al. (2019) investigated the orbital stability and habitability of exomoons around known planets of M dwarfs. They found four planets that were in the habitable zone and for which the orbit of the moon is stable for a long timescale (the migration time of the moon is larger than the Hubble time). These four planets are HIP 12961 b, HIP 57050 b, GJ 876 b and c (in their paper under the names of CD-23 1056b, Ross 1003b, IL Aqr b, and c, respectively). Interestingly, for these planets our model gave 0%, 0%, 6% and 0% survival rates, respectively. One of these four planets, however, HIP 57050 b was thoroughly investigated also by Trifonov et al. (2020) (planet name: GJ 1148b in their paper) where they were using secular theory and direct N-body integrations to study the possibility of a moon. Similar to our results, they have found that it is unlikely that a moon would stay in a stable orbit around the planet.
Zollinger et al. (2017) investigated the tidal evolution and habitability of exomoons assuming possible configurations and found that moons in the habitable zone of M dwarfs experience a tidal heating rate comparable to the tidal heat flux of Jupiter’s moon, Io (~2 W m⁻²). For stars with masses <0.2M☉, because of the strong gravitational effect of the nearby star, the tidal heating rate of a Mars-like moon around a Jupiter-like planet is much stronger than that of Io’s, implying that these moons might not be habitable. The maximum limit of tidal heating for habitability, however, is not known, and atmospheric effects can influence the habitability of bodies, as well (Barnes et al. 2013).

Recent hydrodynamical simulations showed that icy moons can form from circumplanetary disks of ice giants such as Uranus and Neptune, indicating that there might be a large population of moons around exoplanets beyond the snowline (Szulágyi et al. 2018). According to our findings, these moons, if formed, can stay in stable orbits for a long period of time.

Older stars host more evolved systems, in which there is a higher possibility that most satellites have disappeared due to either disruption (at the Roche radius of the planet) or escape. The presence of a massive exomoon may prevent synchronization of the planet with the star and promotes synchronization with the moon. In this case, satellites also evolve into a 1:1 spin-orbit resonance, getting to a mutually synchronised planet-moon system. Less dissipative planets (gas giants and ice giants) are more prone to keep satellites because of slower orbital migration (see e.g., Sasaki et al. 2012). If an exomoon is observed around a planet, it can provide a strong constraint on the value of the quality factor, *Q*, which describes the dissipation in the planet and gives hint on its interior.

Since younger stars host less evolved systems, there is a higher chance for planets to still have satellites (which will escape or be destroyed later). But satellites on less distant orbits are more difficult to observe (Szabó et al. 2006; Simon et al. 2012). If more than one moon is found, then their orbital architecture may give information on their formation, and provide a strong constrain on *Q*. It may also be possible to infer the presence of rings.

## 6. Conclusion

There is no confirmed exomoon detection to date, but the quest for finding moons in Kepler data is still in progress. Teachey et al. (2018) reported an exomoon candidate, Kepler-1625 b I (see also Teachey & Kipping 2018), but this discovery is under debate (Kreidberg et al. 2019; Teachey & Kipping 2019). Regardless of the existence of this moon, we found stable orbits for satellites around this planet. Recently, the presence of an exomoon candidate was suggested based on observations of sodium and potassium lines around planet WASP-49 b (Oza et al. 2019). Further observations are needed for confirmation.

As a result of our calculations, we conclude that planets with long orbital periods support high survival rates for moons, and close-in planets, which are easier to observe, are less likely to host moons. This can serve as an explanation for not having a confirmed exomoon discovery so far.

Thanks to high accuracy photometry and long-term follow-up missions, CHEOPS is measuring the radius of numerous known exoplanets through transits with unprecedented accuracy. This paves the way for characterizing the close environment of exoplanets and can unveil the presence of exomoons and exorings. Moons and rings are natural companions of planets. Neptune-sized planets are especially promising targets owing to their (putative) lower dissipation and large mass compared to Earth-like planets. Such targets will be in the reach of the CHEOPS mission. Following the CHEOPS mission, ESA's PLATO 2.0 mission, capable of discovering and following exoplanet systems and host stars, and the ARIEL mission will continue the task of collecting precise exoplanet light curves. Whereas exomoons have never been firmly detected up to date, they should exist and their successful detections will provide invaluable information on the planet's interior and formation process (Crida & Charnoz 2012). Detecting the first exomoons will be a major discovery as moons could provide habitable environments, in particular, if they orbit volatile-rich planets (like giant planets) and are located inside, or close to, the habitable zone (Forgan & Dobos 2016; Dobos et al. 2017).

We thank the referee, Alessandro Trani for the useful comments which improved the manuscript.

The COFUND project oLife has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 847675. This project has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grant GINOP-2.3.2-15-2016-00003, the City of Szombathely under Agreement No. 67.177-21/2016 and the grant KEP-7/2018 of the Hungarian Academy of Sciences. V.D. has been supported by the Hungarian National Research, Development, and Innovation Office (NKFIH) grants PD-131737 and K-131508. S.C. acknowledges support for the CHEOPS mission from the CNES agency, and EXOATMOS program from LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001).

## ORCID iDs

A. Pál https://orcid.org/0000-0001-5449-2467

## References

- Alvarado-Montes, J. A., Zuluaga, J. I., & Sucerquia, M. 2017, MNRAS, 471, 3019

- Barnes, J. W., & O'Brien, D. P. 2002, ApJ, 575, 1087

- Barnes, R., Mullins, K., Goldblatt, C., et al. 2013, Asbio, 13, 225

- Canup, R. M. 2004, ARA&A, 42, 441
Canup, R. M. 2012, *Sci*, **338**, 1052

Canup, R. M., & Ward, W. R. 2006, *Natur*, **441**, 834

Chambers, J. E., & Wetherill, G. W. 1998, *Icar*, **136**, 304

Charnoz, S., Morbidelli, A., Dones, L., & Salmon, J. 2009, *Icar*, **199**, 413

Chen, J., & Kipping, D. 2017, *ApJ*, **834**, 17

Cowan, N. B., & Agol, E. 2011, *ApJ*, **729**, 54

Crida, A., & Charnoz, S. 2012, *Sci*, **338**, 1196

Dobos, V., Heller, R., & Turner, E. L. 2017, *A&A*, **601**, A91

Domingos, R. C., Winter, O. C., & Yokoyama, T. 2006, *MNRAS*, **373**, 1227

Forgan, D., & Dobos, V. 2016, *MNRAS*, **457**, 1233

Gavrilov, S. V., & Zharkov, V. N. 1977, *Icar*, **32**, 443

Goldreich, P., & Soter, S. 1966, *Icar*, **5**, 375

Guimaraes, A., & Valio, A. 2018, *AJ*, **156**, 50

Heller, R. 2012, *A&A*, **545**, L8

Kipping, D. M. 2009, *MNRAS*, **392**, 181

Kreidberg, L., Luger, R., & Bedell, M. 2019, *ApJL*, **877**, L15

Lainey, V., Arlot, J. E., Karatekin, Ö., & van Hoolst, T. 2009, *Natur*, **459**, 957

Lainey, V., Dehant, V., & Pätzold, M. 2007, *A&A*, **465**, 1075

Lainey, V., Karatekin, Ö., Desmars, J., et al. 2012, *ApJ*, **752**, 14

Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, *Icar*, **281**, 286

MacDonald, G. J. F. 1964, *RvGSP*, **2**, 467

Martínez-Rodríguez, H., Caballero, J. A., Cifuentes, C., Piro, A. L., & Barnes, R. 2019, *ApJ*, **887**, 261

Miguel, Y., & Ida, S. 2016, *Icar*, **266**, 1

Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, *Icar*, **219**, 737

Nagasawa, M., Thommes, E. W., Kenyon, S. J., Bromley, B. C., & Lin, D. N. C. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. of Arizona Press)

Ogilvie, G. I., & Lin, D. N. C. 2004, *ApJ*, **610**, 477

Oza, A. V., Johnson, R. E., Lellouch, E., et al. 2019, *ApJ*, **885**, 168

Rainey, E. S. G., & Aharonson, O. 2006, in 37th Annual Lunar and Planetary Science Conf. 37, ed. S. Mackwell & E. Stansbery (Houston, TX: Lunar and Planetary Institute), 2138

Ray, R. D., Eanes, R. J., & Lemoine, F. G. 2001, *GeoJI*, **144**, 471

Remus, F., Mathis, S., Zahn, J. P., & Lainey, V. 2012, *A&A*, **541**, A165

Sasaki, T., Barnes, J. W., & O'Brien, D. P. 2012, *ApJ*, **754**, 51

Sasaki, T., Stewart, G. R., & Ida, S. 2010, *ApJ*, **714**, 1052

Scholz, A., Moore, K., Jayawardhana, R., et al. 2018, *ApJ*, **859**, 153

Simon, A. E., Szabó, G. M., Kiss, L. L., & Szatmáry, K. 2012, *MNRAS*, **419**, 164

Smith, J. C., & Born, G. H. 1976, *Icar*, **27**, 51

Sucerquia, M., Ramírez, V., Alvarado-Montes, J. A., & Zuluaga, J. I. 2019, *MNRAS*, **489**, 2313

Szabó, G. M., Szatmáry, K., Divéki, Z., & Simon, A. 2006, *A&A*, **450**, 395

Szulágyi, J., Cilibrasi, M., & Mayer, L. 2018, *ApJ*, **868**, L13

Teachey, A., & Kipping, D. M. 2018, *SciA*, **4**, eaav1784

Teachey, A., & Kipping, D. M. 2019, AAS Meeting, 233, 404.03

Teachey, A., Kipping, D. M., & Schmitt, A. R. 2018, *AJ*, **155**, 36

Tokadjian, A., & Piro, A. L. 2020, *AJ*, **160**, 194

Trifonov, T., Lee, M. H., Kürster, M., et al. 2020, *A&A*, **638**, A16

Zhang, K., & Hamilton, D. P. 2008, *Icar*, **193**, 267

Zollinger, R. R., Armstrong, J. C., & Heller, R. 2017, *MNRAS*, **472**, 8