Exomoons in the Habitable Zones of M Dwarfs

As of 2019 March 19, the Extrasolar Planets Encyclopaedia⁶ had reported a total of 4011 exoplanets with the label "confirmed" (Schneider et al. 2011). In our analysis, first, we discarded the 258 planets detected by direct imaging (123), microlensing (88), pulsar timing (33), transit timing variation (8), and astrometry (6), many of which are located at large distances from the solar system and orbital separations, suffer from strong high-energy radiation, are difficult to characterize, or even are bona fide brown dwarfs (Caballero 2018).

Of the remaining 3734 planets, 205 orbit 109 systems with M-dwarf hosts detected via transit and/or radial velocity. To identify them, we selected all planets orbiting stars with M spectral type in the Extrasolar Planets Encyclopaedia except for nine giants and subgiants (HD 208527, HD 220074), pre-main-sequence stars (CVSO 30, V830 Tau), and spectroscopic, eclipsing, and close binaries (DW UMa, HD 41004 B, Kepler-64, KIC 5095269 AB, NLTT 41135). Besides, we selected all planets orbiting stars without spectral type in the Encyclopaedia, but with masses lower than 0.71 M_⊙, which is the largest mass of an M dwarf in the Encyclopaedia, and that have M spectral type in the NASA Exoplanet Archive.⁷ Among the 205 exoplanets, we also included the planets orbiting the trio of stars BD–21 784, BD–06 1339, and Kepler-155 (Lo Curto et al. 2013; Rowe et al. 2014; Tuomi et al. 2014), which have been spectroscopically classified as stars at the K–M boundary (Upgren et al. 1972; Stephenson 1986; Zechmeister et al. 2009; Muirhead et al. 2012a; Gaidos et al. 2013; Tuomi et al. 2014), and LP 834-042 c, which is listed in the Encyclopaedia as detection_type=Other but has a radial velocity confirmation (Astudillo-Defru et al. 2017b).

The 109 M dwarfs and 205 exoplanets are listed in the Appendix tables. For the stars, we tabulate common/discovery and alternative names,⁸ equatorial coordinates and parallactic distances from Gaia DR2 (Gaia Collaboration et al. 2018), except for Luyten's star and Lalande 21185, for which we used Hipparcos parallactic distances (van Leeuwen 2007). Distances range from merely 1.30 pc (Proxima Centauri; Anglada-Escudé et al. 2016) to about 430 pc (Kepler-235; Borucki et al. 2011). All stars further than 30 pc have been discovered by Kepler (and K2) except for NGTS–1 and HATS–6 (Hartman et al. 2015; Bayliss et al. 2018). Besides, most of them are field M dwarfs, but at least four belong to stellar cluster members or associations (namely K2-25 in the Hyades–Mann et al. 2016a, K2-33 in Upper Scorpius–Mann et al. 2016b, K2-95 and K2-264 in the Praesepe–Pepper et al. 2017; Livingston et al. 2019).

For the planets, we tabulate their names, letters, main orbital parameters (semimajor axes a and orbital periods P_orb), and corresponding references. For the 43 planets without a published semimajor axis (e.g., Rowe et al. 2014; Muirhead et al. 2015; Sahlmann et al. 2016), we derived them from their known periods and our estimated stellar masses, using Kepler's Third Law. For the 89 planets discovered with the RV method and without measured transits, we tabulate their minimum masses, ${M}_{2}\sin i$ (marked with ">"), while for the 103 planets discovered with the transit method and without RV follow-up, we tabulate their radii. We have true masses and radii for 13 transiting planets with RV follow-up in 10 systems (see Table 1).

RV measurements only yield the minimum mass of an exoplanet, whereas transit measurements only yield the radius. For our habitability analysis, we used planet masses and radii. However, only 13 planets have both RV and transit measurements, i.e., true masses and radii (see above). As a result, for the other 192 planets, we derived planetary radii and masses. Among the multiple mass–radius relations available in the literature for different planetary compositions (e.g., Fortney et al. 2007; Seager et al. 2007; Sotin et al. 2007; Lissauer et al. 2011; Barnes et al. 2013; Wu & Lithwick 2013; Weiss & Marcy 2014; Barnes et al. 2015; Rogers 2015; Wolfgang et al. 2016; Chen & Kipping 2017; Kanodia et al. 2019), we chose that of Chen & Kipping (2017), as it covers a broad range of planetary sizes.

The largest planet transiting an M dwarf is NGTS–1 b, with R ∼ 1.33 R_Jup (Bayliss et al. 2018). On the other side of the relation, the smallest planets transit the stars Kepler 138 b, LSPM J1928+4437 (b, c, and d), Kepler-125 c, and 2MUCD 12171 (now universally known as TRAPPIST-1; d and h), and have masses and radii in the intervals 0.10–0.39 M_⊕ and 0.52–0.78 R_⊕, close to the values of Mars (0.107 M_⊕, 0.343 R_⊕). However, the distributions of masses and radii of the whole planet sample peak around 5–10 M_⊕ and 1–3 R_⊕, respectively.

First of all, we compiled Gaia DR2 equatorial coordinates (in J2015.5 epoch and J2000.0 equinox) with VizieR (Ochsenbein et al. 2000) for 107 of our stars. For the other two, Lalande 21185 and Luyten's star, without a Gaia DR2 entry, we retrieved the equatorial coordinates from Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), which were separated by less than two months from the J2000 epoch. Next, we collected multiband broadband photometry from the ultraviolet to the mid-infrared from the literature: GALEX DR5 (FUV, NUV; Bianchi et al. 2011), SDSS DR9 (ugri; Ahn et al. 2012), UCAC4 (BV, gri; Zacharias et al. 2013), Pan-STARRS1 DR1 (gri; Tonry et al. 2012; Chambers et al. 2016), Tycho-2 (B_T, V_T; Høg et al. 2000) APASS9 (BV, gri; Henden et al. 2016), 2MASS (${{JHK}}_{s};$ Skrutskie et al. 2006), Gaia DR2 (G_BP G G_RP; Gaia Collaboration et al. 2016, 2018), WISE (W1W2W3W4; Cutri et al. 2012), and allWISE (W1W2W3W4; Cutri et al. 2014).

For that, we used the Upload X-Match tool of the interactive graphical viewer and editor for tabular data TOPCAT (Ochsenbein et al. 2000) with a search radius of 5'' and the "Best" find mode option. We made sure that the automatic cross-match was right by carrying out an inspection of each individual target and datum using the Aladin interactive sky atlas (Bonnarel et al. 2000), which provides simultaneous access to digitized images of the sky, astronomical catalogs, and databases. Eventually, for the 109 stars we compiled 1682 individual magnitudes along with their uncertainties.

We homogeneously derived the bolometric luminosity and effective temperature for each star with the compiled photometric data and the Virtual Observatory Spectral Energy Distribution Analyzer⁹ (VOSA; Bayo et al. 2008). We selected the best fit between a grid of BT-Settl-CIFIST theoretical spectral models (Baraffe et al. 2015) and our empirical data under the constraints log g = 4.5, 2000 K < T_eff < 4500 K, and [Fe/H] = 0.0. The VOSA fits are acceptable in all cases, except for the extremely young star K2-33 in the Upper Scorpius association (David et al. 2016; Mann et al. 2016b; Schweitzer et al. 2019) and the metal-poor star K2-155 (Díez Alonso et al. 2018b; Hirano et al. 2018). Further details on the luminosity calculation will be described by C. Cifuentes et al. (2019, in preparation).

The luminosities of the 109 planet-host stars are displayed in the Appendix. As expected, the most luminous stars are at the K–M boundary, with L ∼ 0.14–0.17 L_⊙. In particular, the most luminous stars, Kepler-155 and Kepler-252, have VOSA-derived T_eff = 4300–4400 K, which are actually warmer than the hottest known M dwarfs (thus, their spectral types should be revised). On the other side, the least luminous stars, with L = (5–15) × 10⁻⁴ L_⊙, are the latest and coolest ones, which in our case are 2MUCD 12171 and Proxima Centauri, which have M8.0 V and M5.5 V spectral types, respectively.

Finally, we computed stellar masses from K_s magnitude and distance using the relations of Mann et al. (2019). We display them in the last column of Table 5. The lowest masses, M_* = 0.09–0.12 M_⊙, are found again for 2MUCD 12171 and Proxima Centauri, and the highest ones, M_* = 0.67–0.70 M_⊙, correspond to the likely late-K dwarfs Kepler-155 and Kepler-252, K2-33 (actually incorrect, due to its extreme youth), and K2-155 (probably due to a metallicity effect). For comparison, we compiled published luminosities for 26 stars from the following works: Leggett et al. (2001), Howard et al. (2010), Boyajian et al. (2012), Gaidos & Mann (2014), Quintana et al. (2014), Hartman et al. (2015), Schlieder et al. (2016), Dittmann et al. (2017), Gillon et al. (2017), Maldonado et al. (2017), Bakos et al. (2018), Hirano et al. (2018), Ribas et al. (2018), Smith et al. (2018), Affer et al. (2019), Díez Alonso et al. (2019), Feinstein et al. (2019), Günther et al. (2019), Hobson et al. (2019), and Perger et al. (2019). The luminosities derived in this work are in good agreement with the values given in the literature (see Figure 1). Only five stars deviate more than 20% of the 1:1 ratio, but only one deviates significantly from our measurement. For BD–17 400, we calculated L = (835 ± 13) × 10⁻⁴ L_⊙, but Leggett et al. (2001) tabulated L = (450 ± 70) × 10⁻⁴ L_⊙ with pre-Gaia data. The other four stars have luminosities consistent within error bars with our measurements (LP 656–38,¹⁰ Leggett et al. 2001; GJ 96, Gaidos & Mann 2014; Kepler-186, Quintana et al. 2014; K2-149, Hirano et al. 2018).

Zoom In Zoom Out Reset image size

We calculated an HZ for each host star using the one-dimensional climate models from Kopparapu et al. (2014). We chose their "recent Venus" and "maximum greenhouse" estimates for the inner and outer HZ, respectively. The exoplanets' oceans might begin to evaporate and condense in the stratosphere at planetary equilibrium temperatures ≈320–340 K ("moist–greenhouse," Kasting et al. 1993; Kopparapu et al. 2013, 2014). However, between this process and the complete oceanic evaporation (runaway greenhouse limit, Kasting et al. 1993; Kopparapu et al. 2013, 2014), the exoplanets could still have surface temperatures adequate for habitability ("continuosly HZ," Hart 1978; Kasting et al. 1993).

Figure 2 illustrates the HZ computation for our compiled exoplanet-host M dwarfs. The exoplanet candidates are depicted in three different colors, depending on whether they lie between their host stars and the inner HZ boundary (146), whether they lie between the inner and outer HZ boundary (33), or whether they lie beyond the outer HZ boundary (26). The 33 potentially habitable exoplanets are marked with the string "Yes" in column "Potentially HZ" from Table 6. Of the 33 planet candidates, 10 have masses or minimum masses above 10 M_⊕ and 9 between 5 and 10 M_⊕. Not counting the inclination angle effect in the RV planets, which translates into larger actual masses, most if not all of these 19 planets are mini-Neptunes and can hardly sustain liquid water on a rocky surface (Fortney et al. 2007; Miller-Ricci et al. 2009; Dorn et al. 2017). In other words, these 19 planets lie within the HZ of their stars, but they are unlikely to be habitable. However, their Earth-sized moons could actually be habitable.

Zoom In Zoom Out Reset image size

As illustrated by Figure 2, three habitable planets are beyond the approximate 1 Gyr-tidal locking boundary: BD–06 1339 c (Lo Curto et al. 2013), Kepler-186 f (Torres et al. 2015), and HD 147379 b (Reiners et al. 2018). Of them, two are Neptune-like planets and one, Kepler-186 f, has a mass of only about 1.8 M_⊕. However, following Barnes (2017), they will likely be tidally locked after 4.6 Gyr, the solar system age (see Section 3.2 for an individualized computation of planet locking times).

We focused on the 33 exoplanets orbiting the HZ of their host stars and analyzed whether or not they could host exomoons on reasonable timescales. For that, we followed the work of Piro (2018), who considered null exomoon eccentricity and obliquity for simplicity. If an exoplanet hosting an exomoon is close enough to its host star, the star's gravity strips the exomoon away. This critical distance, namely the radius of the Hill sphere, is given by:

$$\begin{eqnarray}&&{a}_{\mathrm{crit},m}={f}_{\mathrm{crit}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{* }}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$$

where M_p and M_* are the exoplanet and star masses, and f_crit = 0.49, a critical factor that assumes a prograde orbit (Domingos et al. 2006; Piro 2018). Eventually, the exomoon will be stripped away by the host star if ${a}_{m}\gt {a}_{\mathrm{crit},m}$, or fall back into the planet if ${a}_{m}\lt {a}_{\mathrm{crit},m}$, where a_m is the moon semimajor axis around the exoplanet.

If an exomoon gets too close to its host exoplanet, it will be tidally disrupted upon reaching the Roche limit, i.e., the point where the exomoon can no longer bind its material by gravitational forces. This condition is given by (Piro 2018):

$$\begin{eqnarray}&&{a}_{\mathrm{Roche}}\approx 1.34{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{\rho }_{m}}\right)}^{1/3},\end{eqnarray} \tag{ 2 }$$

where ρ_m is the exomoon's density (see also Aggarwal & Oberbeck 1974). Under the assumption that all the initial angular momentum of an exoplanet eventually goes into its exomoon, the exomoon's semimajor axis can be written as:

$$\begin{eqnarray}&&{a}_{m}=\displaystyle \frac{4{\pi }^{2}{\lambda }^{2}{M}_{{\rm{p}}}{R}_{{\rm{p}}}^{4}}{{{GM}}_{m}^{2}{P}_{0}^{2}},\end{eqnarray} \tag{ 3 }$$

where λ (λ_⊕ = 0.33, Williams 1994) is the exoplanet's radius of gyration (${I}_{{\rm{p}}}/{M}_{{\rm{p}}}{R}_{{\rm{p}}}^{2}$, see Piro 2018), M_m is the exomoon mass, I_p is the planet moment of inertia, R_p is the planet radius, and P₀ is the planet initial spin period. To estimate the exomoon migration timescale, namely, the e-folding timescale for the change in the exoplanet–exomoon distance, we used Equation (25) from Piro (2018), which can be rewritten as:

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig},m}={t}_{k}\left(\displaystyle \frac{{k}_{2,\oplus }}{{k}_{2}}\right)\left(\displaystyle \frac{{\tau }_{\mathrm{lag},\oplus }}{{\tau }_{\mathrm{lag}}}\right)\left(\displaystyle \frac{{a}_{m}^{\,8}}{{M}_{m}\,{R}_{{\rm{p}}}^{\,5}}\right),\end{eqnarray} \tag{ 4 }$$

where t_k = 250 Gyr, and a_m, M_m, and R_p are the exoplanet–exomoon separation, exomoon mass, and planetary radius in units of Earth–Moon separation, Moon mass, and Earth radius, respectively.¹¹ Here, k₂ is a Love number of degree 2 (${k}_{2,\oplus }=0.3$) and τ_lag is the time lag between the passage of the perturber and the tidal bulge (Barnes 2017; Green et al. 2017; Piro 2018). For the Earth, ${\tau }_{\mathrm{lag},\oplus }=638\,{\rm{s}}$ (Lambeck 1977; Neron de Surgy & Laskar 1997). The time lag τ_lag is related with the tidal quality factor Q via (Henning et al. 2009; Matsumura et al. 2010; Barnes et al. 2013):

$$\begin{eqnarray}&&{Q}^{-1}\approx n{\tau }_{\mathrm{lag}},\end{eqnarray} \tag{ 5 }$$

where n is the mean motion, which implies that

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{lag},\oplus }}{{\tau }_{\mathrm{lag}}}=\displaystyle \frac{Q}{{Q}_{\oplus }},\end{eqnarray} \tag{ 6 }$$

where Q_⊕ = 12 (Williams et al. 1978). Therefore, combining Equations (3), (4) and (6),

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig},m}\propto Q\displaystyle \frac{{R}_{{\rm{p}}}^{27}{M}_{{\rm{p}}}^{8}}{{P}_{0}^{16}{M}_{m}^{17}}.\end{eqnarray} \tag{ 7 }$$

Table 2 summarizes the used values of k₂, Q, and λ for four representative planetary types. These parameters are reasonable extrapolations of k₂, Q, and λ from models of exoplanets and actual measurements of planets of different sizes and compositions (Dickey et al. 1994; Yoder 1995; Aksnes & Franklin 2001; Zhang & Hamilton 2006; Henning et al. 2009). For planets with sizes between R_p = 1.50 R_⊕ (small, rocky worlds) and R_p = 3.88 R_⊕ (Neptune), we fitted a linear function to k₂:

$$\begin{eqnarray}&&{k}_{2}={k}_{\mathrm{2,0}}+\gamma \,{R}_{{\rm{p}}}.\end{eqnarray} \tag{ 8 }$$

For Q and λ, we fitted, respectively, two power laws and two linear functions, between R_p = 1.50 R_⊕ and R_p = 3.88 R_⊕, and between R_p = 3.88 R_⊕ and R_p = 10.97 R_⊕ (Jupiter):

$$\begin{eqnarray}&&Q={Q}_{0}\,{R}_{{\rm{p}}}{}^{\delta }\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\lambda ={\lambda }_{0}+\epsilon \,{R}_{{\rm{p}}}.\end{eqnarray} \tag{ 10 }$$

The values resulting from these fittings are listed in Table 3. Next, we computed the parameters a_Roche, a_m, and ${\tau }_{\mathrm{mig},m}$ for the 33 exoplanets in HZ with a hypothetical exomoon. In our calculations, we set M_m = M_Moon and ρ_m = 2500 kg m⁻³, which is a density intermediate between those of Io and Ganymede, but slightly lower than that of the Moon (that has a small iron core formed after the collision between Theia and the proto-Earth). The Roche limit, a_Roche, smoothly depends on ρ_m (Equation (2)). For example, assuming the Earth density ρ_⊕ = 5513 kg m⁻³, a_Roche would only be increased by 23%. Next, we set ${a}_{m}={a}_{\mathrm{crit},m}$ in Equation (4) when the exomoon is stripped (${a}_{m}\gt {a}_{\mathrm{crit},m}$) and used the value given by Equation (3) when it falls back (${a}_{m}\lt {a}_{\mathrm{crit},m}$). We did not calculate ${\tau }_{\mathrm{mig},m}$ if a_m < a_Roche. We remark that, for multiplanetary systems, these estimates may not be entirely accurate because of the impact of the other planets.

Table 3.  Parameters for Deriving k₂ (Equation (8)), Q (Equation (9)), and λ (Equation (10))

Rp	k2,0	γ	Q0	δ	λ0
		$({R}_{\oplus }^{-1})$	$({R}_{\oplus }^{-\delta })$			(${R}_{\oplus }^{-1}$)
≤1.50	0.3	0	102	0	0.37	0
1.50–3.88	−0.456	0.504	14.019	4.846	0.458	−0.059
3.88–10.97	1.5	0	24.600	4.431	0.217	0.003
10.97–26	1.5	0	106	0	0.254	0

Download table as:  ASCIITypeset image

The initial planetary spin P₀ is one of the greatest unknowns in such estimates. It is the orbital period, as opposed to the initial rotation period, which would be set by the last major impact. Figure 3 shows the orbital periods for all the exoplanets in our sample. Minimum and maximum P_orb range from 4.3 hr (K2-137 b; Smith et al. 2018) to 24 yr (BD–05 5715 c; Montet et al. 2014), but they concentrate between 1 and 50 days. The range of P_orb for the exoplanets in HZ is narrower, from 4 to 130 days, with a median value of 31 days. For completeness, the median P_orb for inner and outer HZ planets are 7 days and 502 days, respectively.

Zoom In Zoom Out Reset image size

Figure 4 depicts ${a}_{\mathrm{crit},m}$ versus a_Roche for the 33 potentially habitable planets in our sample. It shows the two most extreme cases: ${P}_{\mathrm{0,1}}={P}_{\mathrm{orb}}$ and P_0,2 = 3.0 hr. On the one hand, P_0,1 = P_orb (${\tau }_{\mathrm{mig},m}\propto Q\ {R}_{{\rm{p}}}^{27}{M}_{{\rm{p}}}^{8}{P}_{0}^{-16}$) is the longest allowed value of P₀ because the exoplanet is tidally locked to its star. It is probably the most realistic case, as the habitable exoplanets in our sample orbit their stars under synchronous rotation. In particular, we computed the tidal locking time t_lock with "constant-phase lag" (CPL, where the phase between the perturber and the tidal bulge is constant and insensitive to orbital and rotational frequencies) and "constant-time lag" (CTL, where the time interval between the perturber's passage and the tidal bulge is constant) models, using the publicly available EqTide¹² code and the Proxima Centauri b system as a template (Barnes et al. 2016; Barnes 2017). We assumed zero eccentricity and obliquity, an initial planetary spin period of 3.0 h, a stellar Love number k_2,* = 1.5, a stellar gyration radius λ_* = 0.1 (λ_⊙ = 0.061, Bouvier 2013), and the planetary parameters k₂, Q, λ varying accordingly with R_p. We chose the minimum t_lock between the calculations from the CPL and the CTL models.

Zoom In Zoom Out Reset image size

We found that only four habitable planets have t_lock between 0.1 Gyr and 1 Gyr, namely LP 834–042 d (Astudillo-Defru et al. 2017b), PM J11293-0127 d (Almenara et al. 2015; Sinukoff et al. 2016), HD 156384 C e (Anglada-Escudé et al. 2013), and CD–44 11909 c (Tuomi et al. 2014), and 9 greater than 1 Gyr, GJ 96 b (Hobson et al. 2018), CD–23 1056 b (Forveille et al. 2011), BD–06 1339 c (Lo Curto et al. 2013), Ross 1003 b (Trifonov et al. 2018), HD 147379 b (Reiners et al. 2018), V1428 Aql b (Kaminski et al. 2018), Kepler-186 f (Torres et al. 2015), and IL Aqr b and c (Trifonov et al. 2018). In this situation, all of them could host exomoons beyond the Roche limit after all the material in the protoplanetary disk is dissipated (10–30 Myr, Bertout 1989; Hartmann et al. 1998; Calvet et al. 2002). All other habitable planets get tidally locked in less than 90 Myr.

However, the moon migration timescales are too low. In this scenario, only for six habitable exoplanets the condition a_m > a_Roche holds. The migration timescales are greater than t_H, being t_H = 14.4 Gyr the Hubble time (which is slightly longer than the universe age at 13.8 Gyr), for four exoplanets (CD–23 1056 b, Forveille et al. 2011; Ross 1003 b, Haghighipour et al. 2010; Trifonov et al. 2018; IL Aqr b and c, Trifonov et al. 2018), whose properties are shown in the top part of Table 4 and discussed in Section 3.3. The four of them have ${a}_{m}\gt {a}_{\mathrm{crit},m}\gt {a}_{\mathrm{Roche}}$, the largest planetary radii in the HZ (R_p = 12–14 R_⊕), and τ_mig ≫ 0.8 Gyr, which was the Earth's Late Heavy Bombardment duration (Fassett & Minton 2013; Norman & Nemchin 2014). In this P_0,1 = P_orb scenario, the hypothetical moons around the other 29 planets would fall onto them on timescales much shorter than protoplanetary disk dissipation times under the EqTide layout, not counting potential viscosity in the protomoon disk (Lin & Papaloizou 1986; Pollack et al. 1996; Montmerle et al. 2006; Sasaki et al. 2010). The fifth and sixth habitable exoplanets for which a_m > a_Roche, LP 834–042 b and BD–06 1339 c, have short τ_mig of about 50 and 18 Myr, respectively, which are comparable to the typical protoplanetary disk dissipation time (Ida & Lin 2004; Hernández et al. 2007; Lambrechts & Johansen 2012).

Table 4.  Potentially Habitable Exoplanets That Can Host Exomoons

Planet	GJ	a	Porb	M*	Mp	Rp	am	aRoche	${a}_{\mathrm{crit},m}$	${\tau }_{\mathrm{mig},m}$	${P}_{\mathrm{orb},m}$	Qtidal
		(au)	(days)	(M⊙)	(M⊕)	(R⊕)	(aM)	(aM)	(aM)	(Gyr)	(days)	(TW)
CD-23 1056 b	...	0.25	53.435	0.62	>114.0	(>13.2)	2.72	0.23	2.72	>tH	≤11.0	≲7
Ross 1003 b	1148	0.17	41.38	0.35	>96.7	(>12.0)	2.08	0.21	2.08	>tH	≤8.0	≲32
IL Aqr b	876	0.21	61.082	0.34	>760.9	(>13.3)	5.36	0.43	5.36	>tH	≤11.8	≲5
IL Aqr c	876	0.13	30.126	0.34	>241.5	(>14.0)	2.29	0.29	2.29	>tH	≤5.9	≲126

Download table as:  ASCIITypeset image

On the other hand, P_0,2 = 3.0 hr (${\tau }_{\mathrm{mig},m}\propto Q\ {R}_{{\rm{p}}}^{27}{M}_{{\rm{p}}}^{8}$) is the smallest value used by Piro (2018), similar to the shortest P_orb in our sample (P_0,2 ≈ 4.13 hr, K2-137 b), which is close to its star and far from the HZ boundary, and probably near the shortest allowed value of P₀ (see below). The value of 3.0 hr is also lower than the sidereal rotation periods of Jupiter and Saturn (9.925 and 10.55 hr, respectively), and might be unrealistic in the context of this work. Contrary to the P_0,1 = P_orb scenario, all moons are stripped away (which means that ${\tau }_{\mathrm{mig},m}\propto Q\ {R}_{{\rm{p}}}^{-5}{M}_{{\rm{p}}}^{8/3}{M}_{* }^{-8/3}$), but with migration timescales greater than 0.8 Gyr for 26 planet–moon systems, including the four systems in Table 4. Of them, 15 planet–moon systems have τ_mig > t_H. They are GJ 96 b (Hobson et al. 2018), CD–23 1056 b (Stassun et al. 2017), LP 834–042 b (Astudillo-Defru et al. 2017b), LP 834–042 d (Astudillo-Defru et al. 2017b), Kapteyn's star b (Anglada-Escudé et al. 2014), BD–06 1339 c (Lo Curto et al. 2013), PM J11293–0127 d (Almenara et al. 2015; Sinukoff et al. 2016), Ross 1003 b (Trifonov et al. 2018), HD 147379 b (Reiners et al. 2018), HD 156384 C e (Anglada-Escudé et al. 2013), CD–44 11909 c (Tuomi et al. 2014), V1428 Aql b (Kaminski et al. 2018), Kepler-186 f (Torres et al. 2015), and IL Aqr b and c (Trifonov et al. 2018). However, except for nine cases, GJ 96 b, CD–23 1056 b, BD–06 1339 c, Ross 1003 b, HD 147379 b, V1428 Aql b, Kepler-186 f, and IL Aqr b and c, all tidal locking times are shorter than 0.8 Gyr (${t}_{\mathrm{lock}}\ll {\tau }_{\mathrm{mig},m}$), which actually transforms this scenario into the previous one, P_0,1 = P_orb. For the four planets listed in Table 4, the migration timescales remain unscathed, as ${a}_{m}\gt {a}_{\mathrm{crit},m}$.

To better understand the values of the exomoon migration timescales, we calculated the maximum and minimum allowed ${\tau }_{\mathrm{mig},m}$, i.e., we substituted ${a}_{m}={a}_{\mathrm{crit},m}$ and a_m = a_Roche in Equation (4), respectively. Figure 5 shows the range of ${\tau }_{\mathrm{mig},m}$ for the potentially habitable exoplanets in our sample. In principle, all the 24 exoplanets above the Earth's Late Heavy Bombardment line might host exomoons. However, the parameter space is only broad enough for the five exoplanets on the upper right side of Figure 5, which are the four planets listed in Table 4 and BD–06 1339 c. We compared our migration timescales for these five most extreme cases with the values resulting from applying the numerical technique of Piro (2018) by varying the planet–moon separation over an allowable range and the spin of the planet, and found that our values are slightly longer than, but comparable to, these. We kept our values for consistency.

Zoom In Zoom Out Reset image size

As presented in Section 3.2, only four potentially habitable exoplanets might be able to host exomoons for timescales above 0.8 Gyr under the P₀ = P_orb scenario: CD–23 1056b (a giant exoplanet around an active, metal-rich, M0.0 V star in Eridanus–Forveille et al. 2011), IL Aqr b,c (two giant exoplanets around a nearby M4.0 V star with a four-planet system—Delfosse et al. 1998; Marcy et al. 1998, 2001; Trifonov et al. 2018), and Ross 1003b (an eccentric, Saturn-mass planet around a nearby M4.0 V star with a multiplanetary system¹³ —Haghighipour et al. 2010; Trifonov et al. 2018). The four planets are relatively large and have minimum masses between 0.30 and 2.39 M_Jup (i.e., between 97 and 761 M_⊕).

We analyzed the effects of tidal heating Q_tidal (e.g., Barnes et al. 2008; Jackson et al. 2008a, 2008b, 2008c) on the four exomoons to assess if their mantles would melt, which, for a rocky body, occurs if Q_tidal ≳ 100 TW (Moore et al. 2007; Lainey et al. 2009; Veeder et al. 2012; Driscoll & Barnes 2015). Barnes & Heller (2013) defined four categories of Earth-mass exoplanets in the HZ depending on the tidal heating that they experience: "Earth twins" (Q_tidal < 20 TW), "Tidal Earths" (20 TW < Q_tidal < 1020 TW), "Super-Ios" (1020 TW < Q_tidal < 134000 TW), and "Tidal Venuses" (Q_tidal > 134000 TW). We estimated Q_tidal for each exoplanet–exomoon system by means of the publicly available VPLanet¹⁴ code (Barnes et al. 2016, 2018, 2019), modifying the example IoHeat. So far, we had assumed zero eccentricity and obliquity, as in Piro (2018), in which case tidal heating can only come from rotation, assuming that it is not synchronous. Hence, we relaxed this condition and set a small value for the eccentricity e = 0.01. We also assumed that the exomoon is tidally locked to the exoplanet. As Q_tidal is directly proportional to the moon's size, we took a Ganymede-like exomoon to estimate an upper limit for the tidal heating (Ganymede is the biggest solar system moon).¹⁵ We set the exomoon's semimajor axis a_m to the mean distance between the Roche limit and the Hill radius from Table 4. We let each system evolve for t = 10⁴ years and found that only a potential exomoon around IL Aqr c might reach the mantle melting threshold. We report the upper limit for ${Q}_{\mathrm{tidal}}$ in Table 4.

Some of these moons might be detected with current and near-future technology. The detectability of moons by the transit method (or, more likely, photocentric transit timing variation) has been extensively and intensively discussed in the literature (Sartoretti & Schneider 1999; Barnes & O'Brien 2002; Szabó et al. 2006; Pont et al. 2007; Kipping 2009; Kipping et al. 2014; Hippke & Angerhausen 2015; Simon et al. 2015; Perryman 2018). The ESA space mission (CHEOPS; Fortier et al. 2014) might discover moons with the transit method around the brightest M dwarfs, while, if the moons were massive and dense enough, ESPRESSO at the Very Large Telescope (Pepe et al. 2010) could also detect the modification of the amplitude and signature of the Rossiter–McLaughlin effect as in Zhuang et al. (2012). Unfortunately, none of the four planets in the stars' HZ with long moon migration times transits. If the NASA space mission TESS, currently in operation (Ricker et al. 2014), discovered additional M-dwarf planets in this class, but transiting, they should be a high priority target for ESPRESSO or any other ultra-stable spectrograph with cm s⁻¹ accuracy at a large telescope.

Regarding the potential habitability of these Moon-mass exomoons, they may not retain an atmosphere at 300 K due to their low gravities. As Titan is less than twice as massive (1.83 M_M), increasing the exomoon mass to Titan's would not significantly affect our results; the condition a_m > a_crit,m for the four systems in Table 4 would still be met. The minimum mass to hold onto an atmosphere in the HZ of a star is not known, but is probably close to Mars' ($8.71\,{M}_{{\rm{M}}}$, Pollack et al. 1987; Selsis et al. 2007). If we increased the masses of our putative exomoons to Mars' mass, then, for the four systems in Table 4, ${\tau }_{\mathrm{mig},m}$ would be of the order of 5 million years for CD-23 1056b and Ross 1003b, and would lie above the Hubble time for IL Aqr b anc c. These four systems would remain near the upper right corner of Figure 5, indicating that their ${\tau }_{\mathrm{mig},m}$ values are sufficiently long to still permit a long lifetime.