Impact of Tides on the Potential for Exoplanets to Host Exomoons

Over 4000 exoplanet candidates have been discovered in the past few decades, yet there has been no confirmed exomoons (although see Teachey & Kipping 2018 for a possible detection). Using Earth's moon as a standard, in our own solar system there are five moons at least as massive, and several within an order of magnitude. The abundance of moons in our neighborhood strongly suggests that exomoons will eventually be found to be ubiquitous in most exoplanet systems.

The presence of a moon is not just a curiosity, since it may play an important role in the habitability of the hosting planet. A moon with mass similar to Earth's moon or larger may give rise to a number of significant effects. It can stabilize the obliquity of a planet so that the axial tilt does not vary erratically (Laskar et al. 1993). This prevents extreme temperature swings at the surface of a planet, which might be crucial for some organisms, depending on their adaptability. Additionally, a moon may break the tidal locking of a planet to its star, thereby preventing a day-only and night-only hemispheric split on the planet. As shown in Piro (2018), a moon also drives tidal heating within a planet. This helps power tectonic activity in the planet crust, which aids in the movement of gases in the carbon cycle and may play a role in the prevention of the runaway greenhouse effect surmised to have occurred on Venus (Jackson et al. 2008). A large enough moon may also prevent the cooling of the planet core, thereby sustaining the dynamo effect that drives the magnetic field (Andrault et al. 2016). The presence of a magnetic field blocks harmful radiation and preserves the atmosphere, in contrast to the fate of Mars. Finally, for those planets that do not have a solid surface, life has the potential to develop on the moon itself.

Given the importance that moons may have, we explore the potential for exoplanets to host exomoons. Using parameterized models for the tidal dissipation in the planet, we follow the spins and orbital separation evolution for planet–moon–star systems. These are applied to known Kepler exoplanets to assess the length of time they could potentially host a moon. This is an extension of previous work by Martínez-Rodríguez et al. (2019) which looked into the dynamical stability of exomoons around M-dwarf exoplanets, but here we compute the full time evolution to understand the moon survival times of single-planet systems. We highlight exoplanets within the habitable zone of their host star which could hold a moon for over approximately a gigayear, the rough timescale for the development of life (Schopf et al. 2018). Many current exoplanet candidates lack a mass measurement, but we demonstrate that the ability for their exoplanets to hold a moon for a significant time depends strongly on the strength of their tidal dissipation. Thus, if these planets were found to host moons, it would provide important information about the planet's structure.

In Section 2, we summarize the planet–moon separations allowed in these systems, and introduce the two tidal lag models used in our analysis. In Section 3, we explore the general properties of these models to show the dependence on different parameters. In Section 4, we apply the models to real exoplanets, and then discuss notable specific systems in Section 5. We describe the consequences and future implications of this study in Section 6 and provide a summary in Section 7.

We first summarize the basic framework we use for studying the tidal interactions for a system of a single planet, a moon, and a star. Such systems have been addressed in a number of previous studies (e.g., Ward & Reid 1973; Touma & Wisdom 1994; Neron de Surgy & Laskar 1997; Barnes & O'brien 2002; Sasaki et al. 2012; Sasaki & Barnes 2014; Adams & Bloch 2016). Here we follow the basic picture presented in Piro (2018). Considering just the planet and the moon, the moon generates two bulges on the planet, one on the close side and one on the far side. However, the planet cannot react instantaneously to the tidal forcing. Thus, the bulge will not directly align with the position of the moon, but will be offset, depending on the composition, rigidity, and dissipation rate of the planet. This is referred to as a tidal lag. If the planet and moon are not tidally locked, the planet will spin faster (or slower) than the moon is orbiting. Thus, the moon will essentially be lagging (or leading) the bulge, causing a net torque that will slow down (or speed up) the rotation of the planet. This in turn will cause the moon to spiral outwards (or inwards).

In the following subsections, we introduce the governing equations that describe this process. We then show the limits on the planet–moon separation parameter and detail the two tidal lag models used in our analysis.

2.1. Governing Equations

To solve for the time evolution of the planet–star–moon systems, we use a set of coupled differential equations describing the rotational frequency Ω_p of the planet, the planet–moon separation a_m, and the planet–star separation a_s. The resulting three differential equations describe the spin evolution of the planet due to the combined torques from the star N_s and the moon N_m,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({I}_{p}{{\rm{\Omega }}}_{p})={N}_{s}+{N}_{m},\end{eqnarray} \tag{ 1 }$$

the orbital evolution of the planet,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}[{M}_{p}{\left({{GM}}_{s}{a}_{s}\right)}^{1/2}]=-{N}_{s},\end{eqnarray} \tag{ 2 }$$

and the orbital evolution of the moon,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}[{M}_{m}{\left({{GM}}_{p}{a}_{m}\right)}^{1/2}+{I}_{m}{n}_{m}]=-{N}_{m}.\end{eqnarray} \tag{ 3 }$$

The masses of the planet, moon, and star are M_p, M_m, and M_s, respectively, and the moment of inertia of planet and moon are I_p and I_m respectively. The moment of inertia of the star is not included because we assume that the inertia is relatively large so that the star's spin is insignificantly altered by the tides of the planet. While this approximation may break down for hot Jupiters whose masses are much closer to the star mass and are in short period orbits, it does not affect the habitable zone planets which are the main focus of this paper. The orbital frequency of the moon is denoted as n_m and the term ${I}_{m}{n}_{m}$ in Equation (3) assumes that the moon is tidally locked to the planet because of the relatively short locking time of the moon, which is on the order of millions to tens of millions of years for the Earth–Moon system. Whether or not the moon is actually tidally locked to the planet has little effect on our main results given the small spin angular momentum of the moon.

This set of equations makes a number of simplifying assumptions. All orbits are taken to be circular, thus we do not include the rate of change of eccentricity. In addition, we consider secular timescales so that the tidal bulges caused by the moon are torqued by only the moon and tidal bulges caused by the star are torqued only by the star and not vice versa. This is because over many orbits the torque from one body on the bulge caused by another will cancel. Although giant planets may undergo significant structural evolution which may impact a moon's tidal evolution at early times (Alvarado-Montes et al. 2017), for this work we assume that the moon survives these initial phases. Finally, we assume that the spin and orbital angular momenta are all aligned.

2.2. Moon Separation Limits

Here we summarize how we set the initial semimajor axis for the moon. The outer limit depends on the distance between planet and star, as the closer in the planet is to its star, the more likely the star is to strip the moon away from the planet. The maximum separation between planet and moon also depends on the ratio between planet mass and star mass, and is given by

$$\begin{eqnarray}&&{a}_{m,\max }={{fa}}_{s}{\left(\displaystyle \frac{{M}_{p}}{3{M}_{s}}\right)}^{1/3}.\end{eqnarray} \tag{ 4 }$$

A constant value of 0.49 is used for f, as detailed in Domingos et al. (2006). For a system in which the moon's orbital rate is less than the planet's spin rate, the moon will slow down the planet's rotation and move further out in its orbit (as is the case currently for Earth). Should the moon exceed this maximum distance, it will be lost to the parent star.

On the other hand, there is also a minimum separation required for stability between planet and moon. This is set by where the planet's tidal forces on the moon exceed the moon's self gravity, the so-called Roche lobe. This inner boundary is thus dependent on the radius of the moon and the ratio between planet mass and moon mass (Frank et al. 2002),

$$\begin{eqnarray}&&{a}_{m,\min }=2.16{R}_{m}{\left(\displaystyle \frac{{M}_{p}}{{M}_{m}}\right)}^{1/3}.\end{eqnarray} \tag{ 5 }$$

In detail, the moon may be subject to tidal disruption or send it colliding with the planet as it reaches this semimajor axis. Indeed, the ultimate fate depends on the relative densities between planet and moon as discussed in more detail in Piro (2018). For this work, we are mainly focused on the lifetime of the moon itself and not on the details of its final fate.

2.3. Tidal Lag Models

As celestial bodies exert tidal forces on one another, their surfaces and interiors are squeezed and pulled, causing deformations that are eventually dissipated. This tidal dissipation is a complex problem explained by a variety of proposals. For instance, Ogilvie (2014) details the process as a fluid dynamics problem caused by a coupling of the gravitational quadrupole moment of one body with the monopole of the other. Other models include dissipative creep processes in rocky material (Storch & Lai 2014) and excitation of oscillatory modes in gas giants (Vick et al. 2019).

Due to the complicated process involved in the dissipation of tides, we consider parameterized models to encapsulate the main impact of tidal interactions. We consider two different formulations for this to better show the range of uncertainty inherent in this approach. Although we find quantitative differences between the different parameterizations, our results are qualitatively the same in each case.

The first is the constant time lag (CTL) model. In this model, we take the delay between gravitational pull and bulge appearance to be a constant time, depending on the structure of the planet. In this case, the torque is proportional to the time lag τ of the planet and the difference in the rotational frequency of the planet and orbital frequency of the planet around the star or moon around the planet. Indeed, since the time it takes for the bulge to appear is constant, the faster the planet rotates relative to how fast the moon orbits, the larger the torque will be. Conversely, as the planet and moon (or planet and star) approach tidal locking the torque shrinks and approaches zero. The equations for the torques detailing this process are (Ogilvie 2014)

$$\begin{eqnarray}&&{N}_{s}=3{k}_{2}\tau ({n}_{s}-{{\rm{\Omega }}}_{p})\displaystyle \frac{{{GM}}_{s}^{2}{R}_{p}^{2}}{{a}_{s}^{6}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{N}_{m}=3{k}_{2}\tau ({n}_{m}-{{\rm{\Omega }}}_{p})\displaystyle \frac{{{GM}}_{m}^{2}{R}_{p}^{2}}{{a}_{m}^{6}}.\end{eqnarray} \tag{ 7 }$$

Here k₂ is the Love number of the planet that depends on its rigidity. We take this k₂ to be 0.3 for rocky planets like Earth (Yoder 1995), 0.38 for gas giants like Jupiter (Gavrilov & Zharkov 1977), and 0.34 for ice giants like Neptune (see the interval described in Kramm et al. 2011). Furthermore, the time lag τ is taken to be 638 s (Lambeck & Runcorn 1977; Neron de Surgy & Laskar 1997) for rocky planets, 0.766 s for Neptune-like planets, and 0.00766 s for Jupiter-like planets. The basis for timescales in the latter two is from taking ratios between the quality factor Q described in the constant phase lag (CPL) model (see below).

An example of applying the CTL model is shown in Figure 1. Here we start with the present Earth, Sun, and Moon system and with these initial conditions, we integrate the coupled differential Equations (1), (2), and (3) forward in time and plot the resulting Earth–Moon separation. For simplicity for illustrative purposes here, we have excluded all other solar system bodies. In this system, Earth is initially rotating faster than the Moon is orbiting. This causes the Moon to torque Earth down and the Moon starts out moving away until eventually it becomes tidally locked to Earth. At this point, the torque due to the Sun becomes dominant, slowing Earth down further. Thus the Moon torques up Earth and spirals inwards. The process continues and leads to a runaway effect of Earth spinning faster and faster and the Moon moving in closer and closer until it is tidally disturbed after about 1250 Gyr. At this point, Earth would have spun up to have a day equal to about 6 hr, which is thought to be early Earth's rotational period.

Zoom In Zoom Out Reset image size

The second tidal lag model we use is the CPL model. In this model, the planet responds to the gravitational perturbation at a constant phase rather than after a constant time. This means that the angle between the resultant bulge and the line joining the centers of the planet and moon (or star) will always be the same. In contrast to the CTL model, the torque in this model is constant, dependent only on the quality factor Q and the sign of the difference between rotational and orbital frequency.

The quality factor is a measure of the efficiency of the tidal dissipation of the planet and is inversely proportional to the angular lag (Efroimsky & Lainey 2007). So, the higher the Q the smaller the angle and less dissipative the planet. In fact, rocky planets like Earth have a Q of about 12, whereas ice giants like Neptune are much higher at 10⁴ and Jupiters higher still at 10⁶ (Goldreich & Soter 1966; Murray & Dermott 1999).

The time lag τ and quality factor Q are related simply by $2{k}_{2}\tau ({n}_{s}-{{\rm{\Omega }}}_{p})={k}_{2}{\sigma }_{s}/Q$, where ${\sigma }_{s}=\mathrm{sgn}({n}_{s}-{{\rm{\Omega }}}_{p})$ (and similarly n_m and σ_m for the Moon). Thus in the CPL model the torque magnitude is dependent only on Q and the direction is given by the rotational and orbital frequencies. The torques are given by (Efroimsky & Makarov 2013)

$$\begin{eqnarray}&&{N}_{s}=\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{2}}{Q}{\sigma }_{s}\displaystyle \frac{{{GM}}_{s}^{2}{R}_{p}^{2}}{{a}_{s}^{6}},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{N}_{m}=\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{2}}{Q}{\sigma }_{m}\displaystyle \frac{{{GM}}_{m}^{2}{R}_{p}^{2}}{{a}_{m}^{6}}.\end{eqnarray} \tag{ 9 }$$

To summarize, the torque is constant in the CPL model and dependent on the difference between the rotational frequency of the planet and the orbital frequency of the moon around the planet in the CTL model. As the moon slows down or speeds up the planet rotation toward synchronous rotation with the moon's orbit, the torque gets weaker and weaker in the CTL model. For this reason, the CTL model tends to produce longer timescales for moon survival. We show the same Earth–Moon–Sun model in Figure 1, this time using the CPL model. As shown, the Moon now lasts a much shorter 68 Gyr before becoming tidally disrupted.

We next apply our model to theoretical planet, star, and moon systems to explore the general trends we expect them to exhibit. We first show that a planet tidally locked to its star can be "unlocked" by the presence of a moon. Consider a single rocky planet with no moon that is tidally locked to its star, that is, its rotational frequency is equal to its orbital frequency around the parent star. We will use the CPL model for this example, although the CTL model gives qualitatively similar results. Figure 2 shows such a two-body system: a planet of mass $M=1.3{M}_{\oplus }$ and radius $R=1{R}_{\oplus }$ with an initial rotational period of 1.4 days orbits a star of mass ${M}_{s}=1{M}_{\odot }$ at an initial distance of 0.85 au and becomes locked to it by approximately 4 Gyr. The plot shows the periods (in days) rather than frequency to be more intuitive. If such a system is left alone, the planet stays tidally locked to its star.

Zoom In Zoom Out Reset image size

Now we add a moon with the mass and radius of our Moon and place it at the midpoint of the minimum and maximum allowed planet–moon separations. We call this separation ${a}_{m,\mathrm{avg}}$. Figure 3 shows the evolution of the now three-body system through time. Due to the tidal effects of the moon, the planet does not lock to its star in this case, but rather temporarily becomes locked to the moon after about 3 Gyr. The planet and moon are not actually locked at this time but follow a sequence of the moon coming in closer and the planet spinning faster until eventually the moon comes close enough to the planet to be disrupted.

Zoom In Zoom Out Reset image size

The climate on a planet is likely very different between a case where it is tidally locked to the star, which would result in one side always being hot and the back side always being cold, and a case where the planet has day/night cycles. This in turn may impact which planets are expected to be favorable for hosting life. In fact, tidally locked planets are more likely to initiate a runaway greenhouse effect or incur atmospheric collapse (Kite et al. 2011), and they are less likely to go through a "snowball planet" phase which is thought to be an important step for the onset of life (Checlair et al. 2017).

There are a plethora of variables that play a role in the tidal evolution models we employ. The parameters of the moon are perhaps the least important. The inner limit of planet–moon separation depends on moon density, which does not vary much when noting the moon densities in our solar system. So for this work we take all moon densities to be $3\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. The torque due to the moon on the planet's bulge depends on moon mass, so to ensure the moon has an influence we take a modest sized moon and set the mass of all moons henceforth to be the same as our Moon's, 1 M_M.

Planet–moon and planet–star separations are crucial parameters, as indicated by the sixth power in Equations (6) and (7). The planet–star separation for exoplanets is relatively well constrained, and the planet–moon separation is a free variable we can choose since we insert an artificial moon in our calculations. Thus, it is important to consider other parameters when choosing an exoplanet list.

To get a sense of which variable is next most important, we make a series of contour plots comparing how long a moon is kept in a three-body system with varying planet mass and radius. We choose the planet to be rocky with density $5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and use the CPL model with Q = 12. The mass of the star and initial planet–star separation are taken to be $1\,{M}_{\odot }$ and 1 au, respectively, for illustrative purposes. We define τ_keep as the timescale the moon is kept by the planet before the moon is lost to the star or falls into the planet. Figures (4), (5), and (6) are contour plots corresponding to τ_keep for different initial rotation period P₀ of the planet. The noise in each of these plots is due to undersampling of the parameter space and is not physical, but the general trends are still evident. In Figure 4, P₀ is 365 days, so the planet is initially locked to its star. In the next two figures we decrease the periods to 1 day (Figure 5) and 12 hours (Figure 6). The results show that the planet radius has a strong impact on the evolution. As we increase radius vertically up on each contour at any constant mass, there is a significant decrease in the time the moon is kept in the system. In contrast, although increasing mass at constant radius increases timescales, the effect is much less significant. For example, in Figure 5 we can take a planet with $M=1{M}_{\oplus }$, $R=1.6{R}_{\oplus }$, which will hold onto its moon for about 8 Gyr according to the contour plot. If we multiply the mass by 1.5 and keep the radius constant ($M=1.5{M}_{\oplus }$, $R=1.6{R}_{\oplus }$), the timescale increases to around 12 Gyr, a 50% increase. Now if we instead multiply the radius by 1.5 and keep the mass constant ($M=1{M}_{\oplus }$, $R=2.4{R}_{\oplus }$), then the timescale drops to under 2.5 Gyr, a near 70% decrease.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The initial rotational period of the planet is important as well, but given the difficulty of constraining exoplanet rotational frequency, we take each planet to be initially locked to its star (following Figure 4) for the rest of this work. Although many known exoplanets are probably not tidally locked to their star, this approximation is a conservative approach since for most cases the moon retention time is actually longer for faster initial rotation rates for the planet. This is because as long as the moon is not lost to the star, there is a longer period of time where the moon is being pushed away, and then the moon has to travel a longer distance coming back in toward the planet. Also, fast rotating Jupiter-like planets have enough angular momentum in principle to push a moon away until it is lost, but given the high mass of such planets, in practice we find that because the tidal forces are relatively weak, it leads to long moon keep times.

We now apply our models to real exoplanet systems filtered from the NASA Exoplanet Archive.³ We choose Kepler exoplanets, but since often the mass is unknown and hence the density as well, the composition of many of these planets is uncertain. For this reason, we create a planet classification system based mainly on radius as shown in Table 1. We first perform a radius cut and split the planets into three groups: rocky Earth-like planets, Neptune-like planets, and Jupiter-like planets. A second density cut is utilized for radii between $1.75{R}_{\oplus }$ and $2.4{R}_{\oplus }$ to differentiate super-Earths and sub-Neptunes. The systems column displays the number of single-planet systems that follow the criteria.

Table 1.  Planet Classification Scheme and Tidal Parameters

Radius $({R}_{\oplus })$	Density $({\rm{g}}\,{\mathrm{cm}}^{-3})$	Type	Systems	Q	k2	$\tau ({\rm{s}})$	Reference
<1.75	Unknown	Rocky	482	12	0.30	638	(1)
1.75–2.4	>3.6	Rocky	331∗	12	0.30	638	(2,3,4,5)
1.75–2.4	2.1–3.6	Rocky or Neptune	329∗	12 or 104	0.30 or 0.34	638 or 0.766	(2,3,4,5)
1.75–2.4	<2.1	Neptune	331∗	104	0.34	0.766	(2,3,4,5)
2.4–3.5	Unknown	Neptune	346	104	0.34	0.766	(1)
3.5–5.7	Unknown	Neptune or Jupiter	99	104 $\mathrm{or}$ 106	0.34 or 0.38	0.766 or 0.00766	(6)
>5.7	Unknown	Jupiter	506	106	0.38	0.00766	(6)

Note. (1) Fulton et al. (2017), (2) Valencia et al. (2007), (3) Fortney et al. (2007), (4) Charbonneau et al. (2009), (5) Odrzywolek & Rafelski (2016), (6) Stern & Levison (2002); ∗from an overlap of 333 unique exoplanets.

Download table as:  ASCIITypeset image

Since for many planets the mass and density are unknown, we allow planets to be treated as two different compositions. For example, the exoplanet Kepler-1638b has a radius of $1.87\,{R}_{\oplus }$, falling under the radius category between $1.75\,{R}_{\oplus }$ and $2.4\,{R}_{\oplus }$, and it has an unknown mass. Thus we consider the mass to be $5.9\,{M}_{\oplus }$ corresponding to a density of $5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and treat it as a rocky planet (super-Earth) and we also consider the mass to be $1.8\,{M}_{\oplus }$ corresponding to a density of $1.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and treat it as a Neptune-like planet (sub-Neptune). (Jupiter-like planets with unknown mass have densities set to $1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$.) These densities were chosen because they are the average densities of the corresponding planet type in our solar system. Allowing for this overlap is important because composition affects moon survivability due to different τ values (for CTL) or different Q values (for CPL). This is significant for two reasons. First, a planet treated as one type may not hold a moon for a relatively long enough time for life to develop, whereas the same planet treated as another type may do so. Thus we increase the likelihood of finding habitable zone planets that meet the criterion of moon survivability over 1 Gyr. Second, if an exomoon is found in the future, we can infer the composition of its host planet by looking at the age of the system and comparing it to the expected moon survivability time. With the overlap considered, we apply our models to a total of 813 rocky planets, 773 Neptune-like planets, and 605 Jupiter-like planets.

Our main result involves running the models on these real exoplanet systems by inserting an Earth's moon-sized moon at ${a}_{m,\mathrm{avg}}=\tfrac{1}{2}({a}_{m,\min }+{a}_{m,\max })$ in each system, integrating the tidal evolution equations forward in time, and then plotting the moon survival time in gigayears against the stellar flux received by the planet in solar units (su). We use the radius, semimajor axis, and when available, the mass of the planet when running the models. When there is no mass or only an upper limit we assign a mass based on the density desired as explained earlier. We use the stellar properties available to calculate the flux received by the planet. Finally, we assume every planet starts out tidally locked to its star, so that its initial rotation period is equal to its orbital period.

Figure 7 shows the first of our results, the CPL model applied to rocky exoplanets. The dark horizontal line is at 1 Gyr, the approximate time required for life to develop on a planet. The two gray vertical lines mark the optimistic habitable zone for G-type stars as described in Kopparapu et al. (2013), ranging from 0.29 to 1.78 su. Thus we are interested in the small rectangle formed by these lines in the upper middle part of the plot. The planets lying in that rectangle are those that are in the habitable zone of their respective systems and can hold a moon for at least 1 Gyr.

Zoom In Zoom Out Reset image size

In Figure 7, we can note that there are three exoplanets in the rectangle we are interested in. These objects are described in more detail in the next section. Figure 8 shows the CTL model applied to the same rocky group of exoplanets. As mentioned earlier, the timescales are longer in this model (for every planet receiving less than 1.78 SU). Indeed we now have seven exoplanets in the rectangle and note that the three planets from the CPL plot have been promoted from blue to green. Past 1.78 su many timescales are actually lowered when moving from CPL to CTL. This is probably because these planets are very close in, so that the initial rotation rate is relatively fast leading to a large torque in the CTL model. Thus although in the CTL model the torque diminishes with time, the large initial torque is enough to compensate so that the overall timescale is larger. However, if we focus on the habitable zone rectangle, we have longer timescales for the CTL model for the planets of most interest.

Zoom In Zoom Out Reset image size

Figures 9 and 10 show the CPL and CTL model, respectively, for Neptune-like exoplanets. Because these types of planets are significantly less dissipative, we have much longer timescales overall. In fact, there are 27 exoplanets in the rectangle of interest, with most of them holding a moon for longer than 15 Gyr. The blue triangles in Figure 8 that are in the rectangle of interest are all green triangles in these plots, showing that the same planet can have very different timescales based on composition. Also, many rocky planets that are in the habitable zone but have a moon survival time of less than 1 Gyr (in between the two vertical lines but under the horizontal line in Figures 7 and 8) are promoted to the rectangle of interest when treated as Neptune-like. Note the shift in the vertical axis for the Neptune-like planet: the minimum time is raised from 10⁻² to 10⁻⁷ yr eliminating the red symbols (moons lost to star) that do exist.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Jupiter-like planets are displayed in Figures 11 and 12. Most of these planets are close in to their parent star, so we see a trend downward in timescales from CPL to CTL. However, because Jupiter-like planets are even less dissipative than Neptune-like planets, the overall timescales are longer. There are 12 planets in the rectangle of interest, all of which have a moon survival time of over 15 Gyr. Another feature of these plots is the significant increase in filled symbols, which are planets that have a mass determination. This is expected since these are large planets whose mass is easier to measure using the radial velocity follow-up measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the next section, we look into more detail at the exoplanets falling into our category of interest.

Our results show that there are at least 36 exoplanets in the habitable zone of their respective systems that can potentially hold a moon for longer than 1 Gyr. These are summarized in Table 2, along with the moon survival time for each model and composition. Also included is the distance the system is from Earth, and the apparent G-band magnitude of the star.

Table 2.  Habitable Zone Planets with ${\tau }_{\mathrm{keep}}$ >1 Gyr

Planet	Flux (SU)	Type	${\tau }_{\mathrm{keep}}$ CPL (Gyr)	${\tau }_{\mathrm{keep}}$ CTL (Gyr)	Distance (pc)	G-band (mag)	TESS Sector
Kepler-452b	1.104	Rocky	8.57	>15	560	13.40	S14,S15
Kepler-1638b	1.382	Rocky or Neptune	1.60, >15	>15, >15	879	14.75	S14,S15
Kepler-22b	1.091	Rocky or Neptune	5.34, >15	>15, >15	190	11.64	S15
Kepler-1544b	0.890	Rocky or Neptune	0.23, 15	2.79, >15	349	14.04	S14,S15
Kepler-1606b	1.402	Rocky or Neptune	0.68, >15	8.55, >15	880	15.11	S14,S26
Kepler-1653b	1.038	Rocky or Neptune	0.05, >15	0.47, >15	755	NA	S14,S15
Kepler-1090b	1.200	Rocky or Neptune	0.91, >15	11.16, >15	879	14.99	S14,S15
Kepler-443b	0.880	Rocky or Neptune	0.29, 9.95	2.66, >15	821	15.90	S15,S26
Kepler-1552b	1.100	Neptune	>15	>15	786	14.93	S14
Kepler-1632b	1.270	Neptune	>15	>15	732	13.07	S14,S15
Kepler-1540b	0.920	Neptune	6.1	>15	247	14.08	S14,S26
Kepler-896b	1.740	Neptune	11.5	>15	896	15.43	S14,S26
Kepler-991b	1.201	Neptune	1.15	6.1	394	15.15	S14,S15
Kepler-1362b	1.181	Neptune	9.75	>15	747	15.82	S14
Kepler-1539b	1.770	Neptune	9.98	>15	765	15.31	S14,S26
Kepler-1545b	1.371	Neptune	>15	>15	745	15.03	S14,S15
Kepler-1058b	1.202	Neptune	4.78	>15	525	15.5	S14,S15
Kepler-1554b	1.170	Neptune	>15	>15	973	15.71	S14
Kepler-1341b	0.985	Neptune	13.43	>15	483	15.02	S14,S26
Kepler-1318b	0.495	Neptune	>15	>15	481	15.24	S14
Kepler-1600b	0.459	Neptune	>15	>15	1060	15.81	S14
Kepler-1593b	1.020	Neptune	>15	>15	650	15.88	S14,S15
Kepler-1634b	0.651	Neptune	>15	>15	628	14.49	S14,S15
Kepler-1636b	0.850	Neptune	>15	>15	2062	15.79	S14,S15
Kepler-1097b	1.141	Neptune	>15	>15	735	15.21	S14,S26
Kepler-1635b	0.401	Neptune or Jupiter	>15, >15	>15, >15	1100	15.69	S14,S26
Kepler-459b	0.400	Neptune or Jupiter	>15, >15	>15, >15	1525	15.30	S14,S15,S26
Kepler-1625b	1.160	Jupiter	>15	>15	2460	15.76	S14,S15
Kepler-453b	0.921	Jupiter	>15	>15	449	13.52	S14,S26
Kepler-1628b	0.560	Jupiter	>15	>15	352	16.95	S14,S26
Kepler-456b	0.581	Jupiter	>15	>15	758	12.69	S14
Kepler-1519b	1.361	Jupiter	>15	>15	1503	14.80	S14,S15
Kepler-34b	1.250	Jupiter	>15	>15	1897	14.83	S14,S15
Kepler-1654b	0.298	Jupiter	>15	>15	578	13.41	S14,S26
Kepler-86b	1.147	Jupiter	>15	>15	347	12.58	S14,S15,S26
Kepler-1647b	0.575	Jupiter	>15	>15	1256	13.58	S14,S15

Download table as:  ASCIITypeset image

Kepler-22b. The closest system with a possible Earth-like planet is Kepler-22, at 587 light years. Kepler-22b has been an object of interest since its discovery in 2011, as it was the first transiting planet found to be in the habitable zone (Borucki et al. 2012). If the planet is of rocky composition, it can hold a moon for just over 5 Gyr using the CPL model and for over 15 Gyr using the CTL model. If it is Neptune-like, it can potentially hold a moon for over 15 Gyr using either tidal lag model. The system is about 4 Gyr old (Safonova et al. 2016), so regardless of the composition or model, it is likely that if Kepler-22b formed an exomoon, it would retain it until the present day. In fact, the Hunt for Exomoons with Kepler (HEK) project has studied Kepler-22b and constrained a satellite mass to less than $0.5{M}_{\oplus }$ with a 95% confidence (Kipping et al. 2013).

Kepler-1638b. Another notable potentially rocky exoplanet is Kepler-1638b, discovered in 2016 (Morton et al. 2016). With a stellar flux of 1.038 su, it is the most similar to Earth in terms of radiation received. However, if it is a rocky planet it can potentially hold a moon for only 51 Myr or 468 Myr in the CPL and CTL models respectively. On the other hand, if it is actually a sub-Neptune, the moon survival time jumps to over 15 Gyr for either model. A relatively old system, it is estimated to be around 7.7 Gyr, so if a moon is found orbiting Kepler-1638b it is a strong indication that the planet is not actually rocky.

Kepler-991b. Kepler-991b, also discovered in 2016 (Morton et al. 2016), is a Neptune-like planet with a large variation in moon survivability times across the two tidal lag models. In the CPL model, this time is 1.15 Gyr whereas in the CTL model it is 6.101 Gyr. The age of the system is estimated to be around 3.47 Gyr. Thus, should a moon be found around Kepler-991b it would provide a better constraint on the tidal dissipation for this planet.

Kepler-1635b. A notable exoplanet we consider as both Neptune-like and Jupiter-like is Kepler-1635b. For both compositions and tidal dissipation models, this planet has the potential to hold a moon for longer than 15 Gyr. This makes it a solid target for follow-up observations for exomoons in the future. In addition, experimenting with moon size and initial rotation rates in our analysis would separate timescales, with Jupiter-like having the longer timescale. Then a good constraint on the age of the system may reveal which type of planet it is. Either way, Kepler-1635b is certainly gaseous, so it is a good candidate for possibly finding life on an exomoon.

Kepler-1625b. An interesting case is the Jupiter-like exoplanet Kepler-1625b, the only exoplanet (at the time of writing) to host an exomoon candidate (Teachey & Kipping 2018). It is in or near the habitable zone and can potentially hold a moon for over 15 Gyr when using a Moon-sized moon. In fact, the detected exomoon is estimated to be Neptune-sized, bringing the system closer to being binary planets as discussed in Cabrera & Schneider (2007). When running similar calculations, we find that a Neptune-sized moon also stays for over 15 Gyr, regardless of the initial rotational period of the planet, showing that such a moon would be stable for a long time. Figure 13 shows the orbital evolution of the system using three different initial planet rotational periods: 287 days (initially locked to star), 1 day, and 12 hr (similar to Jupiter). In each of these cases, the moon orbit has no significant evolution, signifying the long-term retention of such a moon around this planet. Given the stability of the system, this exomoon may even host a submoon, as shown in Kollmeier & Raymond (2018), but also see Rosario-Franco et al. (2020). Thus, although there is much controversy surrounding the existence of this exomoon (see Heller et al. 2019; Kreidberg et al. 2019; and Rodenbeck et al. 2018 for instance), its apparent detection is a step forward into the future of exomoon discoveries.

Zoom In Zoom Out Reset image size

The results of our study can influence future research in several ways. First, we provide numerous targets that are probable places for exomoons to exist. We have pointed out 36 such planets, all of which are in the habitable zone of their respective systems and can hold a moon for longer than 1 Gyr. Thus, as technology advances and observing power increases, observers can refer to our list for potential targets.

Second, because we have allowed for overlap between planet types, our work can help constrain planet composition in the event of a future detection of an exomoon. In fact, the discovery of an exomoon would eliminate at least one of our many parameters, as we have generically added a moon identical to ours in each of these systems. Reworking our analysis would then provide more accurate results.

Finally, this work adds to our understanding of tidal evolution and specifically the two current tidal lag models. Since in many of our cases the CPL and CTL models predict very different outcomes, future observations may support one of these models over the other.

In this paper we considered few of many parameters. As shown in Piro (2018), initial rotational frequency has a significant impact on moon survival time: shortening the spin period of a planet can dramatically increase this timescale. Also, as mentioned earlier, hot Jupiters can affect their host star's spin rate so that the approximation in Equation (2) breaks down. Thus, including different planet spin periods and the star's spin can be explored in future work.

One of the biggest constraints on our models is that we consider only single-planet systems. This ensures we have simple coupled differential equations that give clean results. However, this also cuts out over 2000 exoplanets, many of which are in the habitable zone of their systems. Thus, a future project can be an extension of this work to multiplanet systems. This would take a reconfiguration of our differential equations, or we may need to turn to more advanced programs should that prove too complex.

In addition to number of planets, there are several other factors that can affect the moon retention timescale that we do not consider here. This includes the possible migration of a planet–moon system inwards leading to a smaller ${a}_{m,\max }$ and potential moon loss. The same effect may send the moon into an evection resonance orbit and collision course toward the planet (Spalding et al. 2016). These issues, among other complex effects, should be considered in future work.

Finally, an exciting part of exoplanet research is the observations currently being made by the Transiting Exoplanet Survey Satellite (TESS), which will view 200,000 nearby stars and catalog thousands of planets. The new discoveries that TESS makes will add to the growing list of habitable zone planets and to our list of potential exomoon hosts as well. With ground-based follow-up observations, exoplanet mass will be much better constrained. This has a significant impact on our work since most of the planets we have analyzed do not have a mass estimation. Thus, in the future we will work to include many more planets in our study and incorporate newly discovered constraints, increasing the number of exoplanets potentially holding an exomoon and improving the accuracy of our results.

TESS observations of the 36 cataloged planets, noted by which sector the planet is viewed in, are included in Table 2.

A current project working on detecting exomoons is HEK, which uses transit timing variations (TTVs) and radial velocity methods to search for exomoons (Kipping et al. 2012), and has detected the exomoon candidate around Kepler-1625b. Although there are no confirmed cases, HEK sets the precedence for using TTVs for exomoon detection, and future missions such as ESA's Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) and CHaracterising ExOPlanets Satellite (CHEOPS) will utilize transits to high precision when looking for exoplanets and possibly exomoons (Cessa et al. 2017; Pascale et al. 2018).

A closely related technique is using transit duration variation, which means how long the transit lasts varies from one orbit to another, signaling the presence of an exomoon or another planet (Kipping 2009). This technique can be utilized by TESS, ARIEL, CHEOPS, and the planned PLAnetary Transits and Oscillations of stars (PLATO) satellite mission (Heng et al. 2014).

Additional methods proposed in finding exomoons include direct imaging with spectroastrometry (Agol et al. 2015), which may be possible with the James Webb Space Telescope using instruments specializing in high-contrast imaging (Boccaletti et al. 2015). A novel idea for searching for exomoons is by following Na i and K i signatures, which are signs of a geologically active satellite (Oza et al. 2019). Although the method is mainly for hot Jupiters in ambient plasma and does not apply directly to the habitable zone exoplanets, our study does include Kepler-14b, a hot Jupiter that can potentially hold a moon for over 15 Gyr, and thus may be a good target to apply this approach. Finally, exomoons may also be detected using microlensing. Liebig & Wambsganss (2010) describe the feasibility of such a process using a single-planet–moon–star system, similar to the configuration chosen for our analysis. This method may be possible with PLATO, although it will specialize in transit techniques. Overall, the number of missions and methods aimed at discovering exoplanets and exomoons make the prospects of a confirmed exomoon a near-future possibility (also see the discussion in Sucerquia et al. 2020.)

In this work, we have explored the potential for exoplanets to host exomoons. We first introduced the coupled differential equations that govern the orbital and rotational evolution of planet–star–moon systems. This includes the fact that a moon can only be stable in a limited range of separations from its parent planet. We detailed the two tidal lag models, CPL and CTL, and applied them to an Earth-like system to show the similarity in final result and difference in timescale.

Next, we analyzed theoretical systems. We first showed an example in which a planet initially tidally locked to its star can be unlocked due to the presence of a moon by considering a system with and without a moon. When calculating moon survival time, the separations a_s and a_m are most influential, but we showed that the planet radius is probably the next most significant parameter. We accomplished this by making a series of contour plots that displayed the trends with different parameters of the systems.

We then applied our models to real data by taking single-planet systems from the NASA Exoplanet Archive. We defined our own planet classification scheme and split exoplanets into rocky, Neptune-like, and Jupiter-like. Our main results are Figures 7 through 12, which show ${\tau }_{\mathrm{keep}}$ versus incident flux on the planet, and the area of interest is a moon survival time of over 1 Gyr and stellar flux such that the planet is in the habitable zone. We found a total of 36 unique exoplanets that follow this criteria: 1 solely rocky, 7 rocky and Neptune-like, 18 solely Neptune-like, 2 Neptune-like and Jupiter-like, and 10 solely Jupiter-like. Thus we find that some of the best locations to host life in the galaxy may be the moons of gas giants.

Notable results we discussed include Kepler-22b: a possible rocky planet that is a good target to search for an exomoon; Kepler-1638b: an exoplanet whose composition can be known should an orbiting moon be found; Kepler-991b: a Neptune-like planet that may help measure tidal dissipation rates; Kepler-1635b: a gaseous planet that is an excellent candidate to have a moon with the potential to harbor life; and Kepler-1625b: an exoplanet with a potential exomoon.

In the future, the detection of moons around exoplanets may be an important constraint on exoplanet structure because the tidal dissipation and thus the time to keep a moon can vary greatly for different types of exoplanets. Such constraints can be assisted with more detailed models of the tidal dissipation since in principle the tidal reaction can be different for forcing from the moon and parent star. This should be explored in future theoretical work.

We thank David Kipping and Shreyas Vissapragada for detailed feedback. We also thank Jaime Alvarado-Montes, José Caballero, René Heller, Apurva Oza, Sean Raymond, Jean Schneider, and Johanna Teske for helpful comments. A.T. acknowledges support from the USC-Carnegie fellowship.