Tidally Heated Exomoons around Gas Giants

We consider three different moons with radii equal to $(1,2,4)\cdot {R}_{\mathrm{Io}}$, which represent Io and roughly Mars- and Earth-sized exomoons. The moons are assumed to be spherically symmetric and made of concentric layers with uniform mechanical and thermal properties (Figure 1). Each moon is assumed to have a metallic core of density ${\rho }_{c}$ and radius R_c and a silicate rocky layer with density ${\rho }_{s}$, which is a standard composition for bodies in this size range. The ratio of the surface (R) and core (R_c ) radii is assumed to be the same for the three exomoons and equal to the value estimated for a sulfur-rich Io core, 0.52 (Anderson et al. 1996). We note that the core ratio could be different for other exomoons; however, the core size has a small effect on our conclusions (see Section 2.2).

Zoom In Zoom Out Reset image size

The outer rocky layer consists of an elastic lithosphere where heat is transferred via conduction and heat piping and a convective mantle (see Section 2.3). In the elastic lithosphere, the temperature decreases linearly to the surface temperature T_surf. Inside the convective mantle, the temperature follows an adiabat (Turcotte & Schubert 2014):

$$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dP}}=\displaystyle \frac{\alpha T}{\rho {C}_{c}}.\end{eqnarray} \tag{ 1 }$$

Here α and C_c are the thermal expansivity and heat capacity, respectively (see Table 1). As an example, for a characteristic mantle temperature of 1800 K, the adiabatic temperature gradient is ≈9 K GPa⁻¹.

Table 1. Model Parameters

Parameter	Units	Formula/Value	Definition
Rg	J mol−1 K	8.3144	Universal gas constant
RIo	km	1821	Radius Io
${\rho }_{m}$	kg m−3	3542	Mantle density
${\rho }_{c}$	kg m−3	5150	Core density1
$\tfrac{{R}_{c}}{R}$	⋯	0.5206	Surface core radius ratio1
α	K−1	2·10−5	Thermal expansivity modulus2
k	W m−1 K−1	4	Thermal conductivity2
κ	m3 s−1	10−6	Thermal diffusivity3
Cc	J kg−1 K−1	1142	Heat capacity
Ea	kJ mol−1	300	Activation energy4
${\eta }_{s}$	Pa s	$1\cdot {10}^{16}$	Mantle solidus viscosity
${\mu }_{0}$	G Pa	65	Mantle shear modulus5
${\mu }_{1}$	K	$8.2\cdot {10}^{4}$	Melt shear modulus coefficient 16
${\mu }_{2}$	⋯	40.44	Melt shear modulus coefficient 2
α	⋯	0.3	Andrade coefficient 1
ζ	⋯	1	Andrade coefficient 2
B	⋯	25	Melt fraction coefficient
n	⋯	2–3	Permeability exponent7
γ	⋯	$\tfrac{{b}^{2}{\rm{\Delta }}\rho g}{\tau {\eta }_{l}}$	Melt scale velocity
S♃	W m−2	50.26	Solar irradiation at Jupiter's orbit
AIo	⋯	0.52	Io's bond albedo8

References. (1) Anderson et al. (1996); (2) Sohl & Spohn (1997); (3) Hussmann & Spohn (2004); (4) Karato & Wu (1993); (5) Segatz et al. (1988); (6) Fischer & Spohn (1990); (7) Moore (2003); (8) Simonelli et al. (2001).

Download table as:  ASCIITypeset image

For sufficiently high mantle temperatures, the local mantle temperature can exceed the local solidus temperature at a given depth. When this occurs, a partially molten sublayer is formed (see Figure 1). The average melt fraction ($\overline{\phi }$) of the layer is computed as

$$\begin{eqnarray}&&\overline{\phi }=\displaystyle \frac{\overline{T}-\overline{{T}_{s}}}{\overline{{T}_{l}}-\overline{{T}_{s}}},\end{eqnarray} \tag{ 2 }$$

with $\overline{T}$ being the average temperature of the layer. $\overline{{T}_{s}}$ and $\overline{{T}_{l}}$ are the averages of the pressure-dependent solidus and liquidus temperatures (Takahashi 1990),

$$\begin{eqnarray}&&{T}_{s}=1409.15+134.2P-6.581{P}^{2}+0.1054{P}^{3},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{T}_{l}=2035+57.46P-3.487{P}^{2}+0.00769{P}^{3},\end{eqnarray} \tag{ 4 }$$

with the pressure in gigapascals. When a slope of ${{dT}}_{s,l}/{dP}=10$ K GPa⁻¹ is reached, we assume that the solidus and liquidus temperatures increase linearly with pressure (Reese et al. 1999).

The lithosphere is assumed to behave elastically, and the viscosity of the mantle depends on the temperature following an Arrhenius relation,

$$\begin{eqnarray}&&\eta ={\eta }_{s}\exp \left(\displaystyle \frac{{E}_{a}}{{R}_{g}\overline{{T}_{s}}}\left(\displaystyle \frac{{\overline{T}}_{s}}{\overline{T}}-1\right)\right),\end{eqnarray} \tag{ 5 }$$

where ${\eta }_{s}$ is the viscosity at the solidus temperature, E_a is the activation energy, and R_g is the ideal gas constant. We do not consider the change of activation energy with pressure, which leads to an increase of activation energy with depth (Karato & Wu 1993), and we assume a constant viscosity for the layer. This can result in an underestimation of the activation energy of the mantle. For example, the activation pressure at Mars' mantle varies between 300 KJ mol⁻¹ close to the surface to 540 KJ mol⁻¹ in the mid-mantle (Nimmo & Stevenson 2000). We expect this variation to have a small effect in our first-order modeling.

The presence of melt weakens the mantle and results in a decrease of viscosity. For the sublayer of melt, we parameterize the decrease of viscosity with melt fraction ($\overline{\phi }$) following Moore (2003) and Henning et al. (2009) and the change of shear modulus with melt fraction using the fit of Fischer & Spohn (1990) and Shoji & Kurita (2014) to laboratory experiments made by Berckhemer et al. (1982):

$$\begin{eqnarray}&&\eta ={\eta }_{s}\exp \left(\displaystyle \frac{{E}_{a}}{{R}_{g}\overline{{T}_{s}}}\left(\displaystyle \frac{{\overline{T}}_{s}}{\overline{T}}-1\right)\right)\exp \left(-B\overline{\phi }\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\mu ={10}^{\tfrac{{\mu }_{1}}{\overline{T}-{\overline{T}}_{s}+1600}-{\mu }_{2}}.\end{eqnarray} \tag{ 7 }$$

The value of B can range from 10 to 40 for a strong or weak mantle, respectively (Henning et al. 2009). We adopt an intermediate value B = 25 consistent with laboratory experiments (Mei et al. 2002). For the layers without melt, we assume a constant shear modulus ${\mu }_{0}$. ${\mu }_{1}$ is an empirical constant, and ${\mu }_{2}$ is adjusted so that the shear modulus is continuous at the solidus temperature.

The melt fraction increases with temperature until the disaggregation point is reached. When this occurs, the asthenosphere does not behave as a viscoelastic solid anymore, and it should be modeled as a magma ocean. We assume that this occurs when $\overline{\phi }\gt 0.45$ (e.g., Moore 2003). Although we do not model heat transport and tidal dissipation in this regime, we discuss this aspect further in Sections 2.2 and 2.3.

We consider two mechanisms of internal heat generation: radiogenic heating and tidal dissipation. Radiogenic heating is computed assuming chondritic composition of the mantle (Schubert et al. 1986) and that the age of the body is the same as the age of the solar system.

Due to tidal interactions with the planet, regular moons evolve into a 1:1 spin–orbit resonance (become tidally locked) soon after formation (e.g., Peale 1999). If the orbit is eccentric or the moon has a nonzero obliquity, tidal forces result in periodic deformations of the moon. As the moon is not perfectly elastic, energy is dissipated in the interior. The total amount of tidal dissipation is (e.g., Cassen et al. 1980; Segatz et al. 1988; Makarov & Efroimsky 2014)

$$\begin{eqnarray}&&\dot{E}=-{\mathfrak{I}}({k}_{2})\displaystyle \frac{{\left({nR}\right)}^{5}}{G}\left(\displaystyle \frac{21}{2}{e}^{2}+\displaystyle \frac{3}{2}{\sin }^{2}\theta \right),\end{eqnarray} \tag{ 8 }$$

where R is the radius of the body, n is the orbital frequency, G is the universal gravitational constant, e and θ are the moons' eccentricity and obliquity, and ${\mathfrak{I}}({k}_{2})$ is the imaginary component of the k₂ Love number, which is related to the quality factor (Q) as $| {\mathfrak{I}}({k}_{2})| =| {k}_{2}| /Q$. We assume that the moon is in a Cassini state with small obliquity, as is the case for Io (Baland et al. 2012), and neglect the contribution of obliquity to the thermal budget. We note that the previous equation is accurate to second order in eccentricity. For high orbital eccentricities, terms of order $O({e}^{4})$ should be added. High-order terms start to have a noticeable effect for $e\gt 0.1$, and orders-of-magnitude differences can arise for very high eccentricities ($e\gt 0.6$; Renaud et al. 2021). Additionally, for $e\gt 0.1$, higher-order spin–orbit resonances might occur (Makarov 2012; Walterová & Běhounková 2020; Renaud et al. 2021), which results in nonsynchronous rotation and increased tidal dissipation before the moon's orbit is circularized.

The value of ${\mathfrak{I}}({k}_{2})$ depends on the internal structure of the body and the tidal frequency. It can be computed using the viscoelastic theory for self-gravitating bodies (Peltier 1974; Wu & Peltier 1982; Sabadini et al. 2016). By using the correspondence principle (Peltier 1974), the equations of motion governing the deformation of each layer can be transformed to the Fourier domain and written as a set of differential equations of the form

$$\begin{eqnarray}&&\displaystyle \frac{d\tilde{{\boldsymbol{Y}}}}{{dr}}=A\tilde{{\boldsymbol{Y}}}.\end{eqnarray} \tag{ 9 }$$

$\tilde{{\boldsymbol{Y}}}$ is a vector containing the Fourier-transformed radial and tangential displacements (${\tilde{y}}_{1}$, ${\tilde{y}}_{2}$), the radial and shear stress (${\tilde{y}}_{3}$, ${\tilde{y}}_{4}$), the gravitational potential (${\tilde{y}}_{5}$), and the potential stress (${\tilde{y}}_{6}$), and A is a matrix given in Appendix A. We use the matrix propagator method of Jara-Orué & Vermeersen (2011) and Sabadini et al. (2016) to solve Equation (9) under appropriate boundary conditions (Appendix A).

To use the correspondence principle, a rheological law relating stress and strain is needed. Different rheological models have been developed for the study of tidally active bodies (e.g., Renaud & Henning 2018); the viscoelastic behavior of the material depends on its shear modulus μ and viscosity η. The simplest and most commonly used is the Maxwell model. The response is characterized by the so-called Maxwell time ($\eta /\mu$). When the tidal period is close to the Maxwell time, tidal dissipation is enhanced. For much shorter forcing periods, the body behaves as an elastic body, while for much longer forcing periods, the body responds as viscous fluid. While the Maxwell model has been widely used for the study of tidally active bodies (e.g., Segatz et al. 1988; Fischer & Spohn 1990; Moore 2003; Henning et al. 2009), it does not properly capture the complex behavior of olivine observed in laboratory experiments (Jackson & Faul 2010; McCarthy & Castillo-Rogez 2013). In particular, the Maxwell model does not incorporate the anelastic transient creep deformation mechanism that describes the viscoelastic behavior of the material over timescales shorter than the Maxwell time, which are relevant for tidal dissipation. The Andrade rheology model (Andrade & Trouton 1910) has been particularly successful in capturing this behavior and adopted in recent studies of tidally active bodies of the solar system (e.g., Castillo-Rogez et al. 2011; Bierson & Nimmo 2016) and other planetary systems (e.g., Walterová & Běhounková 2017; Renaud & Henning 2018). These studies have shown the dramatic influence that changing the rheology from Maxwell to Andrade can have on tidal dissipation models. We use the more realistic Andrade rheology as a baseline for this study and briefly compare it with the classic Maxwell model to illustrate the differences between the two and the implications they have for the thermal–orbital evolution of a moon (Section 2.4.1 and Appendix D).

The Fourier-transformed shear modulus can be defined in terms of the creep function $\tilde{J}$:

$$\begin{eqnarray}&&\tilde{\mu }={\tilde{J}}^{-1}.\end{eqnarray} \tag{ 10 }$$

For the Andrade rheology, the Fourier transform of the creep function ($\tilde{J}$) is (Efroimsky 2012)

$$\begin{eqnarray}&&\tilde{J}=\displaystyle \frac{1}{\mu }-\displaystyle \frac{i}{\eta n}+\displaystyle \frac{{\mu }^{\alpha -1}}{{\left(i\zeta \eta n\right)}^{\alpha }}\alpha !.\end{eqnarray} \tag{ 11 }$$

The first two terms of Equation (11) correspond to the elastic and steady-state creep response characteristic of the Maxwell model, and the last term accounts for the transient creep response included in the Andrade model. Here ζ and α are two empirical parameters that characterize the transient creep response. The value of α for olivine is constrained to vary from 0.1 to 0.5 (e.g., Gribb & Cooper 1998; Jackson et al. 2004). We use an intermediate value of 0.3. The value of ζ depends on the ratio between the anelastic and Maxwell timescales. It can have a dependence on forcing frequency, as the dominant anelastic deformation mechanism might change depending on the forcing period. However, accurate characterization of the deformation mechanisms in this regime is not available, and we make the assumption that diffusion dominates the anelastic response (e.g., Shoji & Kurita 2014; Renaud & Henning 2018; Walterová & Běhounková 2020). In that case, both timescales are equal and $\zeta \approx 1$ (Efroimsky 2012). We caution that for high frequencies or high-stress situations, other deformation mechanisms, such as boundary sliding and dislocation unpinning, might become dominant, resulting in an underestimation of tidal dissipation.

The value of ${\mathfrak{I}}({k}_{2})$ depends on the core size, which can differ for different exomoons. However, its effect is small. For an Io-sized exomoon, dissipation can vary by a factor of 2 for a core size range of 0.1–0.5 if a constant-viscosity mantle is assumed. If a low-viscosity asthenosphere is present, the core size has a smaller effect and changes the dissipation by less than 5%.

When the disaggregation point is reached, the viscoelastic theory for tides cannot be applied. Previous studies of tidal dissipation (e.g Fischer & Spohn 1990; Moore 2003; Renaud & Henning 2018) assumed that when the disaggregation point is reached, viscosity is reduced to that of the liquidus, and dissipation drops dramatically. However, dynamic tides in magma oceans can produce dissipation rates equal to or higher than those observed on Io (Tyler et al. 2015). Tidal dissipation in a magma ocean depends on weakly constrained parameters, such as the magma-ocean dissipation timescale. Instead of explicitly modeling tidal dissipation in the magma ocean, we assume that once the disaggregation point is attained, tidal dissipation remains constant.

The heat generated in the interior of the body via radiogenic heating and tidal dissipation is transported through the mantle of the moon to the surface and radiated to space. While most studies on tidally heated moons and planets consider heat transport via convection only (e.g., Fischer & Spohn 1990; Hussmann & Spohn 2004; Henning et al. 2009; Renaud & Henning 2018), heat transport via heat piping has been proposed to be the heat transport mechanism prevalent in Io (O'Reilly & Davies 1981; Moore 2001, 2003) and Earth's early history (Moore & Webb 2013). Here we consider both convection (Section 2.3.1) and heat piping (Section 2.3.2) as possible heat transfer mechanisms. We do not model heat transfer in the magma-ocean regime but briefly discuss it in Section 2.3.3.

2.3.1. Mantle Convection

Due to the strong dependence of viscosity on temperature, most terrestrial planets in the solar system are in the stagnant-lid convection regime. The only notable exception is Earth, where plate tectonics provide a more efficient way to remove internal heat (Schubert et al. 2001; Korenaga 2013). To estimate the heat transport due to mantle convection, we assume that the exomoon is in the stagnant-lid regime. Tidal heating might trigger plate tectonics (Zanazzi & Triaud 2019). In such a case, the cold near-surface material is recycled, making heat transfer between 20 and 100 times more efficient than heat transport in the stagnant-lid regime and resulting in lower mantle temperatures (Nimmo & Stevenson 2000).

In the stagnant-lid regime, the body's interior is divided into a well-mixed adiabatic convective layer beneath a much stiffer conductive layer (Schubert et al. 1979; Morris & Canright 1984; Solomatov 1995; Reese et al. 1999). The conductive layer is subdivided into a thermal boundary layer, where viscosity decreases by around 3 orders of magnitude, and an immobile stagnant lid. We assume the conductive layer (thermal boundary layer+lid) to behave elastically (see Section 2.1). The amount of heat transported via convection in the mantle can be parameterized as

$$\begin{eqnarray}&&q=k\displaystyle \frac{{\rm{\Delta }}T}{{D}_{m}}{Nu},\end{eqnarray} \tag{ 12 }$$

where k is the thermal conductivity of the mantle, ${\rm{\Delta }}T$ is the temperature increase within the convective layer, D_m is the thickness of the convective layer, and Nu is the Nusselt number. In the stagnant-lid regime, the Nusselt number depends on the Rayleigh number (Ra) and the Frank–Kamenetskii parameter θ:

$$\begin{eqnarray}&&{Ra}=\displaystyle \frac{g\alpha \rho {{D}_{m}}^{3}{\rm{\Delta }}T}{\kappa \eta ({T}_{m})},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\theta =\displaystyle \frac{{\rm{\Delta }}{{TE}}_{a}}{{R}_{g}}.\end{eqnarray} \tag{ 14 }$$

For Newtonian viscosity, scaling arguments lead to (Solomatov 1995; Reese et al. 1999)

$$\begin{eqnarray}&&{Nu}\approx a{\theta }^{-4/3}{{Ra}}^{1/3},\end{eqnarray} \tag{ 15 }$$

with $a\approx 0.5$ being a nondimensional parameter (Reese et al. 1999). Inserting Equations (14)–(15) into Equation (12), the heat flux transported via convection can be obtained (Nimmo & Stevenson 2000; Shoji & Kurita 2014):

$$\begin{eqnarray}&&{q}_{m}\approx \displaystyle \frac{k}{2}{\left(\displaystyle \frac{\rho g\alpha }{\kappa \eta }\right)}^{1/3}{\gamma }_{c}^{-4/3},\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}&&{\gamma }_{c}=\displaystyle \frac{{E}_{a}}{{R}_{g}{T}_{m}^{2}}.\end{eqnarray} \tag{ 17 }$$

For typical mantle conditions, ${\gamma }_{c}$ is of the order of $\approx 0.01$ (Nimmo & Stevenson 2000). The thickness of the elastic lithosphere (D_l ) is computed as

$$\begin{eqnarray}&&{D}_{l}=k\displaystyle \frac{{T}_{m}-{T}_{\mathrm{surf}}}{{q}_{m}},\end{eqnarray} \tag{ 18 }$$

where T_surf is the surface temperature. In the heat-pipe regime, the thickness of the lithosphere and its thermal structure differs from that of the stagnant-lid regime (Kankanamge & Moore 2019; Spencer et al. 2020); the effect on the total tidal dissipation is expected to be negligible.

2.3.2. Heat Piping

As the mantle temperature increases and melt starts to form (Section 2.1), heat piping becomes a more efficient heat transport mechanism than convection. Melt is segregated from the solid mantle and advected upward due to its positive buoyancy. An equilibrium can be reached in which the amount of melt advected upward is compensated for by the melt generated by tidal heating. Moore (2001) developed a model to study heat piping in Io that can be used to compute how much heat can be transported via heat piping inside a rocky body. We follow the approach of Moore (2001) and Bierson & Nimmo (2016) and compute heat transport via heat piping by solving mass conservation and using Darcy's law for porous media (see Appendix B) in the sublayer of melt. The amount of heat transport via heat piping depends on two parameters: the permeability exponent n and the scale velocity γ. Here γ depends on the grain size (b), the density contrast between melt and solid matrix (${\rm{\Delta }}\rho$), the surface gravity (g), the melt's viscosity (${\eta }_{l}$), and a constant (τ) closely linked to the permeability exponent:

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{{b}^{2}{\rm{\Delta }}\rho g}{\tau {\eta }_{l}}.\end{eqnarray} \tag{ 19 }$$

Heat transport efficiency increases with increasing γ and decreasing n. The value of n is typically taken to be between 2 and 3 (Katz 2008). For Io, Moore (2001) estimated γ to range between 10^–5 and 10^–6. We consider both a low- and a high-efficiency melt transport scenario with n = 3 and $\gamma /g\,={10}^{-6}/{g}_{\mathrm{Io}}$ and n = 2 and $\gamma /g={10}^{-5}/{g}_{\mathrm{Io}}$, respectively. By doing so, we are assuming that grain sizes, melt viscosity, and density contrasts are similar for the range of body sizes considered in this study.

2.3.3. Convection in a Magma Ocean

2.4.1. Stable and Unstable Equilibrium States

The thermal state of an exomoon depends on the balance between heat generated in the interior of the body and removed through the lithosphere. The evolution of the mantle temperature can be modeled via a simple equation:

$$\begin{eqnarray}&&{V}_{m}\rho {C}_{c}\left(1+{St}\right)\displaystyle \frac{{dT}}{{dt}}={\dot{Q}}_{\mathrm{int}}-{\dot{Q}}_{s}.\end{eqnarray} \tag{ 20 }$$

V_m is the mantle's volume, St is the Stefan number, ${\dot{Q}}_{s}$ is the surface heat flow, and ${\dot{Q}}_{{\rm{int}}}$ is the internal heat production which includes radiogenic and tidal heating. The surface heat flow is the sum of the heat flow transported via convection and heat piping; coupling between the two mechanisms is neglected. We neglect other energy sources, such as primordial heat and heat due to the cooling of the core. Neglecting the feedback between interior and orbital evolution for now, we can say that when the total heat production equals the total heat lost, an equilibrium state is reached. A thermal equilibrium state is stable if, for a deviation in mantle temperature, the system tends to restore its equilibrium: $d{\dot{Q}}_{s}/{dT}\gt d{\dot{Q}}_{{int}}/{dT}$. If the amount of internally generated heat stays constant, the object can remain in this equilibrium state for a long period of time. The number of equilibrium points depends on the rheological model, as well as the prevalent heat transport mechanism (e.g., Ojakangas & Stevenson 1986; Renaud & Henning 2018). As mentioned before, we use the Andrade model as a baseline and compare it with the classic Maxwell model to underline the differences with previous studies of thermal–orbital evolution of rocky moons (e.g., Fischer & Spohn 1990; Hussmann & Spohn 2004; Moore 2003).

Figure 2 illustrates the location and stability of equilibrium points for an Io-like exomoon for two different mantle rheology (Maxwell and Andrade) and heat transport (convection and heat piping) mechanisms. For the Maxwell model, tidal dissipation is negligible for low mantle temperatures; as the mantle temperature increases and approaches the melting temperature, tidal dissipation sharply increases. A further increase of mantle temperature leads to the formation of a sublayer of melt. Viscosity further decreases within this layer and therein causes a further increase of tidal dissipation. Tidal dissipation peaks at T_res when the Maxwell time of the asthenosphere equals the forcing period. Finally, additional warming of the mantle induces a decrease of viscosity and shear modulus and detunes the asthenosphere from the forcing period. For the Andrade rheology, a similar behavior is observed at high mantle temperatures, but for low mantle temperatures, the transient creep mechanism results in higher heat generation.

Zoom In Zoom Out Reset image size

Equilibrium points for both models are also indicated in Figure 2. In both cases, there is a stable equilibrium point. Moreover, the Maxwell plus convection model has an additional unstable point at $T\lt {T}_{\mathrm{res}}$. A change in the moon's orbital parameters that reduces the amount of tidal dissipation has very different consequences depending on whether the system is at a stable or an unstable point. In the first case, the moon evolves smoothly from one equilibrium point to the other; in the second, the moon enters a runaway cooling phase. As mentioned earlier, the previous analysis holds if the evolution of the orbit is not considered. When the feedback between interior properties and orbital dynamics is considered, a more complex picture arises. As we will show in Appendix D, the differences in the location and stability of the equilibrium points discussed above for the two models have important consequences for the orbital evolution of a moon/exomoon.

In the following subsections, we adopt the more realistic Andrade model and apply it first to Io and then to exomoons to obtain the range of orbital parameters for which hot states can be attained.

2.4.2. The Case of Io

Before applying our model to the more general case of an exomoon, we explore its implications for Io and assess whether it can successfully explain its thermal state.

Figure 3(a) shows Io's heat flux and tidal dissipation as a function of mantle temperature. Depending on melt fraction and the prevalent heat transport mechanism, we distinguish three regimes: the stagnant-lid regime, which corresponds to a mantle without a sublayer of melt; the heat-pipe regime, where melt advection dominates heat transfer; and the magma-ocean regime. Tidal dissipation is computed at Io's current orbital period and different eccentricities. When the eccentricity is low, radiogenic heat is the only significant heat source, and thermal equilibrium is reached for a low mantle temperature, $\approx 1200$ K. As eccentricity increases, the equilibrium state moves to higher mantle temperatures. For moderate eccentricities ($\sim {e}_{\mathrm{Io}}/2$), the dislocation creep mechanism can lead to situations where $d{\dot{Q}}_{s}/{dT}\approx d{\dot{Q}}_{{int}}/{dT}$. In such cases, a quasi-equilibrium state can be attained in which temperature varies slowly and the thermal evolution of the system stagnates (Renaud & Henning 2018). For higher eccentricities, new equilibrium configurations appear at higher mantle temperatures. The mantle temperature and stability of such points depend on the considered heat transfer mechanism (see Figure 3).

Zoom In Zoom Out Reset image size

If we consider Io's current eccentricity, we find that heat transported via mantle convection in the stagnant-lid regime is insufficient to remove Io's present heat flux, as also found in Moore (2003). In contrast, the equilibrium points obtained for the heat-pipe regime are much closer to Io's observed heat flux. For the high-efficient heat transport scenario, we obtain that an equilibrium point is attained for a model with a 40 km thick asthenosphere with a melt fraction of 0.02. The equilibrium heat flux in this configuration is roughly half of Io's. In the low-efficiency heat-piping model, the equilibrium point is reached at roughly two times Io's observed heat flux for a model with a 330 km asthenosphere with a 0.2 melt fraction. The truth probably lies in between, but the second scenario is more consistent with the interpretation of Galileo's magnetometer observations as evidence of a near-surface partially molten layer with a 0.2 melt fraction (Khurana et al. 2011). However, this claim has recently been challenged (Roth et al. 2017; Blöcker et al. 2018). For this model, tidal dissipation is mainly focused in the asthenosphere ($\approx 90 \%$) and is higher in equatorial regions (Figure 4). Assuming that all of the molten rock travels to the surface through channels in the lithosphere, we can estimate the resurfacing rate, which is on the order of 1 cm yr⁻¹. This rate is higher than the minimum 0.02 cm yr⁻¹ required by the lack of impact craters and consistent with the 0.4–14 cm yr⁻¹ resurfacing rate estimated from surface changes observed during Galileo's mission (Phillips 2000).

Zoom In Zoom Out Reset image size

2.4.3. Volcanic Exomoons and Tidally Boosted Surface Temperatures

We apply the approach used to obtain thermal equilibrium states for Io to the more general case of exomoons of different sizes and orbital parameters. We are primarily interested in estimating the surface heat flux and temperature. The average surface temperature can be estimated by considering that the moon is a blackbody with surface temperature T_surf,

$$\begin{eqnarray}&&{T}_{\mathrm{surf}}^{4}=\displaystyle \frac{1}{\sigma }\left(\displaystyle \frac{(1-A)S}{4}+{q}_{s}\right),\end{eqnarray} \tag{ 21 }$$

where σ is the Stefan–Boltzmann constant, S is the stellar irradiation, and A is the moon's bond albedo. As we are considering volcanic rocky worlds around cold exoplanets, we use Io's albedo A_Io and the solar flux at Jupiter's orbit S_♃. Localized volcanic activity can result in spatial and temporal variations of surface temperature; our estimations of surface temperature should be understood as the total thermal output of a rocky exomoon but keeping in mind that how and when this energy is released will depend on the thermal regime of the moon. Moreover, the outgassing of material and formation of a substantial atmosphere would alter the heat balance of the moon and thus the surface temperature (e.g., Noack & Rivoldini 2017).

We can also estimate the amount of outgassed material ($\dot{{M}_{o}}$) and resulting column density. This quantity is closely related to the surface heat flux. A proxy for the outgassing rate can be obtained by considering Io's outgassing rate, ${\dot{M}}_{o,\mathrm{Io}}$, and a tidal efficiency, ${\eta }_{T}$,

$$\begin{eqnarray}&&{\eta }_{T}={\left(\displaystyle \frac{R}{{R}_{\mathrm{Io}}}\right)}^{2}\displaystyle \frac{q}{{q}_{\mathrm{Io}}},\end{eqnarray} \tag{ 22 }$$

as ${\dot{M}}_{o}\sim {\eta }_{T}{\dot{M}}_{o,\mathrm{Io}}$ (Oza et al. 2019; Quick et al. 2020). The previous expression gives an order-of-magnitude estimate; other factors, such as the style of heat transport and volcanism and the mantle composition, are likely to affect the outgassing rate.

We start by considering exomoons of sizes $2{R}_{\mathrm{Io}}$ and $4{R}_{\mathrm{Io}}$ orbiting at Io's orbital frequency (Figures 3(b) and (c)). Some differences with Io (Figure 3(a)) are apparent. For bigger exomoons, higher values of tidal dissipation are attained for the same orbital eccentricity. This should not come as a surprise, as tidal dissipation has a strong dependence on body size ($\dot{E}\propto {R}^{5}$). In contrast, the heat transported via mantle convection has a smaller dependence on body size (Equation (12)). Heat transport via melt advection quickly becomes the dominant heat transport mechanism as orbital eccentricity increases.

If we consider Io's orbital eccentricity, we find surface heat fluxes more than 1 and 2 orders of magnitude higher than Io for Mars- and Earth-sized exomoons (i.e., $2{R}_{\mathrm{Io}}$ and $4{R}_{\mathrm{Io}}$), respectively. By using the previous proxies, we can infer that these exomoons would likely feature more volcanic activity than Io and produce a plasma torus with column densities 2–3 orders of magnitude higher than those observed in the Jovian system for a $2{R}_{\mathrm{Io}}$ and $4{R}_{\mathrm{Io}}$ moon, respectively. Additionally, for $R=4{R}_{\mathrm{Io}}$, we observe that thermal equilibrium is reached for a surface heat flux above solar irradiation. In such circumstances, the second term in Equation (21) becomes important, and tides have a signature in the surface temperature that increases to ∼300 K.

Figure 5 shows the surface heat flux and temperature for different combinations of orbital eccentricity and distance. The orbital distance is shown as a function of the Roche limit, $a/{a}_{R}={\left(2\pi {\rho }_{m}G\right)}^{1/3}/{\left(3{n}^{2}\right)}^{1/3}$, which does not depend on the mass of the orbiting planet. The equilibrium states are obtained considering the high-efficiency heat-piping transport scenario. If the low-efficiency scenario is considered, similar results follow, but the magma-ocean regime is attained for lower eccentricities and orbital frequencies. As expected, surface heat flux and temperature increase with increasing eccentricity and decreasing orbital distance (Equation (12)). For instance, Io would have a surface temperature of 400 K provided it orbited with its current eccentricity at $a/{a}_{R}\approx 2$ or at its current orbital distance with an eccentricity of $\approx 0.1$. The subspace of orbital parameters for which surface heat fluxes are higher than those of Io and tides have an effect on surface temperatures reflects that tidal dissipation increases with body size.

Zoom In Zoom Out Reset image size

From Figure 5, it is clear that, provided an exomoon has a high enough eccentricity and/or orbits close to the planet, super-Ios with intense volcanism and even tidally boosted surface temperatures arise. A question immediately follows: do we expect exomoons to orbit within these regions, and for how long?

As we are considering cold gas giants that, in contrast to close-in exoplanets, are far from their star, we neglect the effect of stellar tides in both the planet and the moon. In an isolated moon–planet system, tidal dissipation within the moon removes energy from the moon orbit and circularizes it. Additionally, tidal dissipation within the planet produces a phase lag in the planets' response to the tide raised by the moon. As a consequence, the moon exerts a torque on the planet that changes the planet's spin rate, and the planet exerts a torque of the same magnitude to the moon that drives orbital migration. The direction of moon migration depends on whether the orbital period of the moon is higher or lower than the rotational period of the planet. We consider that the planet's rotation period is lower than the moon's orbital period. This is justified by the fast rotation rates of the solar system gas giants and those measured for extrasolar gas giants (Snellen et al. 2014; Bryan et al. 2018). Under the previous assumptions, the change of orbital frequency (n), orbital eccentricity (e), and the planet's spin rate ($\dot{\theta }$) are given by (Appendix C)

$$\begin{eqnarray}&&\displaystyle \frac{1}{n}\displaystyle \frac{{dn}}{{dt}}=\displaystyle \frac{1}{{\tau }_{n}}+\displaystyle \frac{{p}_{e}{e}^{2}}{{\tau }_{e}},\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{e}\displaystyle \frac{{de}}{{dt}}=-\displaystyle \frac{1}{{\tau }_{e}},\end{eqnarray} \tag{ 23b }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\dot{\theta }}{{dt}}=-\displaystyle \frac{3}{2}\displaystyle \frac{{{GM}}_{m}^{2}{R}_{p}^{3}}{\kappa {M}_{p}{a}^{6}}\displaystyle \frac{| {k}_{2,p}| }{{Q}_{p}},\end{eqnarray} \tag{ 23c }$$

with ${p}_{e}=57/7$ and where ${\tau }_{n}$ and ${\tau }_{e}$ are the orbital migration and circularization timescales and are given by

$$\begin{eqnarray}&&\begin{array}{l}{\tau }_{n}^{-1}=-\displaystyle \frac{9}{2}\displaystyle \frac{{M}_{m}}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{5}| {\mathfrak{I}}({k}_{2,p})| n\\ \,=\,-\displaystyle \frac{9}{2}\displaystyle \frac{{M}_{m}}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{5}\displaystyle \frac{| {k}_{2,p}| }{{Q}_{p}}n,\end{array}\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&\begin{array}{l}{\tau }_{e}^{-1}=\displaystyle \frac{21}{2}{\left(\displaystyle \frac{{R}_{m}}{a}\right)}^{5}\displaystyle \frac{{M}_{p}}{{M}_{m}}| {\mathfrak{I}}({k}_{2,m})| n\\ \,=\,\displaystyle \frac{21}{2}{\left(\displaystyle \frac{{R}_{m}}{a}\right)}^{5}\displaystyle \frac{{M}_{p}}{{M}_{m}}\displaystyle \frac{| {k}_{2,m}| }{{Q}_{m}}n.\end{array}\end{eqnarray} \tag{ 24b }$$

m and p stand for moon and planet, respectively.

Tidal dissipation within the moon results in orbit circularization and inward orbital migration. The rate at which this occurs depends on the internal structure of the moon via the second-degree Love number ${\mathfrak{I}}({k}_{2,m})$. On the other hand, under the aforementioned assumptions, tidal dissipation within the planet results in outward orbital expansion and despinning of the planet. Combining Equations (23a) and (23c), the change in planet spin rate produced as the moon migrates can be obtained. Unless the moon is very massive and formed very close to the planet, the effects of the moon on the planet's rotation are small. A Mars-like exomoon around Jupiter would slow down Jupiter less than 1% as it migrates from 2a_R to 20a_R . In what follows, we assume $d\dot{\theta }/{dt}\approx 0$, the planet is treated as an infinite source of energy that drives orbit migration. The migration rate depends on how energy is dissipated within the planet, as given by the imaginary component of the planet's second-degree Love number, ${\mathfrak{I}}({k}_{2,p})$. In classical tidal theory, it is considered to be independent of the forcing frequency; we make this assumption and consider a range of frequency-independent ${\mathfrak{I}}({k}_{2,p})$. However, new observations indicate a strong dependence of ${\mathfrak{I}}({k}_{2,p})$ on forcing frequency (Lainey et al. 2017, 2020), which has been linked with the excitation of internal waves in the planet's gaseous envelope (e.g., Ogilvie & Lin 2004). The ${\mathfrak{I}}({k}_{2,p})$ spectrum depends on weakly constrained parameters such as the structure and composition of the planet, which in turn can undergo significant changes as the planet evolves. Implicitly, our model assumes that the exomoon does not excite any of the resonant modes in the planet's dissipation spectrum. The incorporation of such effects would require the analysis of a broad range of properties of the giant planet's dissipation spectrum, as no robust bottom-up model of this spectrum can be set up for poorly constrained exoplanet interiors.

Figure 6 shows the circularization timescales as a function of orbital distance for a hot exomoon with a ${\mathfrak{I}}({k}_{2,m})$ similar to that of Io (Lainey et al. 2009). ${\tau }_{e}$ strongly depends on orbital distance (${\tau }_{e}\propto {\left(a/{a}_{R}\right)}^{13/2}$); if the moon is far from the planet ($a/{a}_{R}\gtrapprox 20$), the circularization timescale is more than 1 Gyr. At such orbital distances, tidal dissipation is quite low. However, a highly eccentric ($e\gt 0.1$) Earth-like exomoon would still experience high levels of tidal dissipation that could even increase the surface temperature to around 500 K (Figure 5). Can we imagine a plausible scenario where this could occur? A possible candidate is the capture of a terrestrial-sized planet by a gas giant via the binary-exchange capture process (Williams 2013). The capture process results in highly eccentric orbits compatible with high values of tidal dissipation.

Zoom In Zoom Out Reset image size

For close-in exomoons, the eccentricities required for high tidal activity are lower (Figure 5). However, the circularization timescales are much shorter, less than a million years for a moon orbiting closer to the planet than Io (Figure 5). A sporadic boost in eccentricity could result in high values of tidal dissipation and boost the surface temperature, but such a boost would inevitably be short-lived, giving little chance for observation. For these moons to exhibit vigorous geological activity, it is necessary that the eccentricity be continuously forced. This can occur via MMRs.

MMRs occur when two moons consistently apply a periodic gravitational perturbation to each other. This happens when the orbital frequencies of the two objects (n₁, n₂) are related via $(p+q){n}_{2}-{{pn}}_{1}-q{\dot{\omega }}_{1}\approx 0$, where ${\omega }_{1}$ is the longitude of pericenter, and p and q are two integers. The periodic gravitational forcing alters the orbit of both objects and excites the eccentricity of the moons. As already mentioned, we cannot explore all possible orbital resonances. Thus, we focus on the $2:1$ orbital resonance ($p=1,q=1$).

In a 2:1 orbital resonance, the resonant variable $v={n}_{1}-2{n}_{2}$ is close to zero. Using the perturbing potential up to first order in eccentricity, the equations governing the orbital evolution of the two moons can be obtained (Yoder & Peale 1981):

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{1}=\displaystyle \frac{{n}_{1}}{{\tau }_{n1}}+\displaystyle \frac{{e}_{1}^{2}{n}_{1}}{{\tau }_{e1}}\left((3+{p}_{e})+3\displaystyle \frac{{\tau }_{e1}}{{\tau }_{e2}}\displaystyle \frac{{M}_{1}}{{M}_{2}}{\left(\displaystyle \frac{{C}_{2}}{{C}_{1}}\right)}^{2}{\alpha }^{1/2}\right)\\ \,+\,3\displaystyle \frac{{M}_{2}}{{M}_{p}}\alpha {\left(\displaystyle \frac{{n}_{1}}{v}\right)}^{2}\dot{v}\left({C}_{1}{e}_{1}-{C}_{2}{e}_{2}\right),\end{array}\end{eqnarray} \tag{ 25a }$$

$$\begin{eqnarray}&&\begin{array}{l}(1+K({e}_{1},{e}_{2}))\dot{v}={n}_{1}\left[\displaystyle \frac{1}{{\tau }_{n1}}\left(1-\displaystyle \frac{{\tau }_{n1}}{{\tau }_{n2}}\right)\right.\\ \,+\,\displaystyle \frac{{e}_{1}^{2}}{{\tau }_{e1}}\left((3+{p}_{e})+12\displaystyle \frac{{M}_{1}}{{M}_{2}}{\alpha }^{2}+\displaystyle \frac{{\tau }_{e1}}{{\tau }_{e2}}{\left(\displaystyle \frac{{M}_{1}}{{M}_{2}}\right)}^{2}{\left(\displaystyle \frac{{C}_{2}}{{C}_{1}}\right)}^{2}\right.\\ \,\times \left.\left.\left(2(6-{p}_{e}){\alpha }^{5/2}+3\displaystyle \frac{{M}_{2}}{{M}_{1}}{\alpha }^{1/2}\right)\right)\right],\end{array}\end{eqnarray} \tag{ 25b }$$

$$\begin{eqnarray}&&{e}_{1}=-\displaystyle \frac{{M}_{2}}{{M}_{p}}\displaystyle \frac{{n}_{1}}{v}\alpha {C}_{1},\quad {e}_{2}=\displaystyle \frac{{M}_{1}}{{M}_{p}}\displaystyle \frac{{n}_{1}}{v}{\alpha }^{3/2}{C}_{2}.\end{eqnarray} \tag{ 25c }$$

Here C₁ and C₂ are two constants with values of −1.19 and 0.428, respectively; α is the ratio between the inner and outer moon semimajor axes, which, close to the 2:1 resonance, equals 0.63; and $K({e}_{1},{e}_{2})$ is a positive number,

$$\begin{eqnarray}&&\begin{array}{l}K({e}_{1},{e}_{2})=-6\displaystyle \frac{{M}_{p}}{{M}_{2}}\displaystyle \frac{1}{{C}_{1}}\left[{\alpha }^{1/2}+{\left(\displaystyle \frac{{C}_{2}}{{C}_{1}}\right)}^{2}\displaystyle \frac{{M}_{1}}{{M}_{2}}\alpha \right]{e}_{1}^{3}\\ \,+\,12\displaystyle \frac{{M}_{p}}{{M}_{1}}\displaystyle \frac{1}{{C}_{2}}\left[1+{\left(\displaystyle \frac{{C}_{1}}{{C}_{2}}\right)}^{2}\displaystyle \frac{{M}_{2}}{{M}_{1}}{\alpha }^{-1/2}\right]{e}_{2}^{3},\end{array}\end{eqnarray} \tag{ 26 }$$

which is small provided the eccentricities are also small.

Different orbital migration timescales (${\tau }_{n1}\ne {\tau }_{n2}$) can lead to the convergence of a pair of moons and the assembly of an MMR (Equation (25a), Figure 7). This can occur either in the protoplanetary disk via differential orbital decay (Figure 7(a)) or, once the protoplanetary disk has dissipated, via differential tidal expansion of the orbit (Figures 7(b) and (c)). In the second case, ${\tau }_{n}$ and ${\tau }_{e}$ are given by Equation (24). In the first case, the orbital evolution is driven by the interactions of the protoplanetary disk with the recently formed moons (Peale & Lee 2002; Canup & Ward 2002). The satellites excite density waves in the protoplanetary disk that cause a torque on the moons and result in inward migration. The protoplanetary disk also dampens the satellite eccentricity. Equation (23) can still be used to model the effects of moon protoplanetary disk interactions in the orbital evolution, provided we use p_e = 3 and define ${\tau }_{n}$ and ${\tau }_{e}$ as (Goldreich & Schlichting 2014)

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{n}}\sim \left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right)\left(\displaystyle \frac{{M}_{d}}{{M}_{p}}\right){\left(\displaystyle \frac{a}{h}\right)}^{2}n,\end{eqnarray} \tag{ 27a }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{e}}\sim \left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right)\left(\displaystyle \frac{{M}_{d}}{{M}_{p}}\right){\left(\displaystyle \frac{a}{h}\right)}^{4}n,\end{eqnarray} \tag{ 27b }$$

where h/a is the aspect ratio of the disk, and M_d is the disk mass.

Zoom In Zoom Out Reset image size

As the moons move deeper into the orbital resonance, the orbital eccentricity increases (Equation (25b)). Depending on the ratio of the migration timescales of the two moons (${\tau }_{n1}/{\tau }_{n2}$), two scenarios are possible: the moons (1) cross the orbital resonance or (2) get trapped in it. The first scenario results in a boost of orbital eccentricity that is later circularized in a timescale ${\tau }_{e}$ (Figure 7(a)). In the second scenario (Figures 7(b) and (c)), an equilibrium eccentricity is reached that can persist for a long time.

For MMR assembled via tidal expansion of the orbit once the protoplanetary disk is dissipated, the first scenario can occur provided the outer moon migrates faster than the inner moon. Given the strong dependence of ${\tau }_{n}$ with orbital distance (Equation (24b)), this requires that $| {\mathfrak{I}}({k}_{2,p}({n}_{1}))| \gt | {\mathfrak{I}}({k}_{2,p}({n}_{2}))|$, which can occur under the resonance-locking scenario (e.g., Fuller et al. 2016; Lainey et al. 2020). The boost in orbital eccentricity due to the crossing of the $2:1$ MMR can be estimated as (Dermott et al. 1988)

$$\begin{eqnarray}&&{e}_{1}={\left(\displaystyle \frac{2\sqrt{6}\tfrac{{M}_{2}}{{M}_{p}}\alpha | {C}_{1}| }{1+4\tfrac{{M}_{1}}{{M}_{2}}{\alpha }^{2}}\right)}^{1/3}.\end{eqnarray} \tag{ 28 }$$

We note that similar expressions can be obtained for different values of p and q; thus, Equation (28) gives an order-of-magnitude estimation of the orbital eccentricity boost for other MMRs. Using the previous equation, we can obtain the boost in orbital eccentricity for different values of ${M}_{1}/{M}_{2}$ and ${M}_{1}/{M}_{p}$. As shown in Figure 8, the boost of orbital eccentricity due to the crossing of a resonance can be close to 0.1. For example, for ${M}_{1}/{M}_{2}=1$ and ${M}_{1}/{M}_{p}={10}^{-4}$, $e\approx 0.05$. Such an increase in orbital eccentricity could result in vigorous and tidally boosted surface temperatures if resonance crossing occurs when the moon is close to the planet. As an example, a Mars-sized exomoon orbiting between Io and Europa's orbit would experience surface temperatures up to 400 K (Figure 5), but the eccentricity would be dampened in less than 10 Myr (Figure 6).

Zoom In Zoom Out Reset image size

As mentioned before, in the second scenario, the two moons remain caught in an orbital resonance for a long period of time. As the moons move deeper into the orbital resonance, the orbital eccentricity increases (Equation (25a)). The increase of orbital eccentricity results in a phase lag between both moons, which prompts a transfer of angular momentum from one moon to another. As we will see, under certain conditions, this transfer of angular momentum ensures that both moons migrate at the same rate ($\dot{v}=0$), and an equilibrium eccentricity is reached. The value of the equilibrium eccentricity can be computed as (Equation (29))

$$\begin{eqnarray}&&\begin{array}{l}{e}_{1}^{2}=-\displaystyle \frac{{\tau }_{e1}}{{\tau }_{n1}}\left(1-\displaystyle \frac{{\tau }_{n1}}{{\tau }_{n2}}\right)\left((3+P)+12\displaystyle \frac{{M}_{1}}{{M}_{2}}{\alpha }^{2}+\displaystyle \frac{{\tau }_{e1}}{{\tau }_{e2}}{\left(\displaystyle \frac{{M}_{1}}{{M}_{2}}\right)}^{2}\right.\\ \,\times {\left.{\left(\displaystyle \frac{{C}_{2}}{{C}_{1}}\right)}^{2}\left(2(6-P){\alpha }^{5/2}+3\displaystyle \frac{{M}_{2}}{{M}_{1}}{\alpha }^{-1/2}\right)\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 29 }$$

The value of the equilibrium eccentricity is proportional to $\sqrt{\tfrac{{\tau }_{e}}{{\tau }_{n}}}$, while the prefactor multiplying this ratio differs depending on p and q. The proportionality of the equilibrium eccentricity to this ratio remains for different values of p and q (Dermott et al. 1988, Equation (39)).

For the assembly of a $2:1$ MMR in the protoplanetary disk, the existence of such equilibrium eccentricity requires that ${M}_{2}\gt 2{M}_{1}$. The value of the equilibrium eccentricity is $e\sim h/a$. For a disk aspect ratio of the order of ∼0.1 (Peale & Lee 2002; Canup & Ward 2002), orbital eccentricities of ∼0.1 can be attained (Figure 9). Once the protoplanetary disk dissipates, the eccentricity is dampened on a timescale ${\tau }_{e}$, or the system relaxes to the equilibrium configuration of the tidally driven scenario.

Zoom In Zoom Out Reset image size

In case orbital migration is driven by tidal forces, the existence of an equilibrium eccentricity requires that the inner moon migrates faster than the outer moon, ${\tau }_{n1}\lt {\tau }_{n2}$, which, for a frequency-independent ${\mathfrak{I}}({k}_{2,p})$, implies ${M}_{2}\lt 20{M}_{1}$. The forced eccentricity is (Equations (24) and (29))

$$\begin{eqnarray}&&{e}^{2}\sim {\left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right)}^{2}{\left(\displaystyle \frac{{R}_{p}}{{R}_{m}}\right)}^{5}\displaystyle \frac{{\mathfrak{I}}({k}_{2,p})}{{\mathfrak{I}}({k}_{2,m})}.\end{eqnarray} \tag{ 30 }$$

Figure 9(a) shows the forced eccentricity for an Io-like moon in a 2:1 orbital resonance with a second moon. We consider that tidal dissipation is negligible in the outer moon and use ${\mathfrak{I}}({k}_{2,p})$ and ${\mathfrak{I}}({k}_{2,m})$ compatible with those estimated for Jupiter and Io (Lainey et al. 2009). For the planet radius, we use the empirical mass–radius relation of Bashi et al. (2017). If the inner-to-outer moon mass ratio is too low, the outer moon migrates too fast, and the moons do not get caught in the MMR. On the other hand, if the ratio is high, the forced eccentricity of the inner moon is small. For a mass ratio equal to that of Io and Europa, we obtain a forced eccentricity of approximately half of Io's present eccentricity. This discrepancy is because Io is part of a more complex resonance chain, the Laplace resonance, which can excite higher eccentricities (Yoder & Peale 1981). By considering the simple 2:1 MMR, we obtain an order-of-magnitude estimation of the forced eccentricity that can be attained if the moon is part of a more complex resonance chain. We further consider the case where both moons have equal mass and obtain the forced eccentricity for different moon and planet masses and different values of ${\mathfrak{I}}({k}_{2,p})$ using Equation (29) (Figure 9(b)). Moderate-to-high orbital eccentricities (${10}^{-3}-{10}^{-1}$) that result in high surface heat fluxes are attained, provided that (1) the moon-to-planet mass ratio is high and/or (2) the planet is highly dissipative. Configurations with ${M}_{m}/{M}_{p}\gt {10}^{-4}$ should be regarded with caution, as satellite formation models for giant planets predict ${M}_{m}/{M}_{p}\sim {10}^{-4}$ (Canup & Ward 2006; Heller & Pudritz 2015a, 2015b). This makes Mars-sized exomoons plausible around super-Jovian (1–12 M_♃) planets, but the ratio for Earth-sized is only attained well inside the brown dwarf regime (>12 M_♃).

The value of ${\mathfrak{I}}({k}_{2,p})$ is uncertain and a matter of intense research. Equilibrium tidal theory suggests $| {\mathfrak{I}}({k}_{2,p})| \sim {10}^{-13}$ (Goldreich & Nicholson 1977); however, astrometric observations of the Jovian and Saturnian systems indicate $| {\mathfrak{I}}({k}_{2,\unicode{x02643}})| \,\sim {10}^{-5}$ and as high as $| {\mathfrak{I}}({k}_{2,\unicode{x02644}})| \sim {10}^{-3}$. High values of $| {\mathfrak{I}}({k}_{2,p})|$ lead to a high equilibrium eccentricity. However, high $| {\mathfrak{I}}({k}_{2,p})|$ implies fast orbital migration (Equation (24)). As the moon migrates outward, tidal dissipation rapidly decreases (Equation (8)), limiting the longevity of a super-Io.

This becomes evident if we study the coupled thermal–orbital evolution of a close-in exomoon. As the formation of a pair of Earth-sized exomoons seems unlikely from moon formation theory, and high values of tidal dissipation are more easily attained for Mars-sized than Io-sized exomoons, we focus on the case of a pair of close-in ($a=2{a}_{R}$) Mars-sized exomoons orbiting a super-Jovian planet (5 M_♃). As shown before, the orbital and thermal evolution of the moon are coupled via the imaginary component of the second-degree Love number (${\mathfrak{I}}({k}_{2,m})$), which in turn depends on the interior structure. We compute tidal dissipation and heat transfer as explained in Section 2.1 and evolve the interior using Equation (20). The orbital evolution is computed using Equation (25). The system is integrated forward in time using a two-step Euler method with an adaptive step size. We assume $| {\mathfrak{I}}({k}_{2,p})|$ is frequency-independent, and we vary its value from 10⁻³ to 10⁻⁶ and study the thermal–orbital evolution of the exomoons. In all cases, we assume that the moons start deep into the resonance with an eccentricity of 0.01.

The evolution of the exomoons for different values of $| {\mathfrak{I}}({k}_{2,p})|$ is depicted in Figure 10. At the start of the simulation, the high eccentricity combined with the close proximity of the exomoon to the planet results in high tidal dissipation and tidally boosted surface temperatures. At this time, the moon is in the magma-ocean regime. The moon's eccentricity changes until the equilibrium eccentricity is reached. The value of the equilibrium eccentricity depends on the value of $| {\mathfrak{I}}({k}_{2,p})|$, with higher values leading to higher values of orbital eccentricity. Meanwhile, the moon migrates at a rate dependent on ${\mathfrak{I}}({k}_{2,p})$, tidal dissipation decreases, and the moon cools down enough to enter the heat-pipe regime. Upon cooling, the value of ${\mathfrak{I}}({k}_{2,m})$ changes, and, in a range of 1–10 Myr, a new equilibrium orbital eccentricity is attained.

Zoom In Zoom Out Reset image size

The close proximity of the moon to the planet combined with high orbital eccentricity results in high surface heat fluxes, as well as tidally boosted surface temperatures. For instance, for $| {\mathfrak{I}}({k}_{2,p})| ={10}^{-4}$, surface heat fluxes of more than 1000 W m⁻² are attained. Surface temperatures are higher than 400 K for the first 10 Myr. This phase is, however, short-lived due to the fast orbital migration of the moon. Tidal dissipation quickly decreases, the melt fraction in the asthenosphere diminishes, and the moon shifts from the magma-ocean regime to the heat-pipe regime. After 100 Myr, surface temperatures are down to 200 K, and after 500 Myr, the contribution of tides to surface temperature is negligible. As the semimajor axis increases, the migration timescale (${\tau }_{n}$; Equation (24)) decreases, but the eccentricity remains nearly constant. A surface heat flux 10–100 times higher than Io is maintained during the first billion years and then decreases to values similar to Io while the moon stays in the heat-pipe regime.

If the planet is less dissipative, the equilibrium eccentricity is lower, but the orbital migration timescale is reduced. These two factors partly compensate for each other; while the eccentricity attained for $| {\mathfrak{I}}({k}_{2,p})| ={10}^{-5,-6}$ is lower than that for $| {\mathfrak{I}}({k}_{2,p})| ={10}^{-3,-4}$, the moon stays closer to the planet for a longer time. As in the more dissipative cases, high surface temperatures are attained during the first million years, but the contribution of tides to surface temperature is negligible after 500 Myr. The surface heat flux stays above that of Io for more than 2 billion yr, giving ample time for the outgassing of material and the formation of a secondary atmosphere and plasma torus.

While we observe small-amplitude fast oscillations of the eccentricity at the beginning of our simulation, our models do not feature the pronounced periodic oscillations characteristic of the models of Hussmann & Spohn (2004) and Fischer & Spohn (1990) for Io. As we show in Appendix D, this is a consequence of the very efficient transfer mechanism included in our model (heat piping) and the use of Andrade rheology, which implies that Io could be in thermal equilibrium instead of an oscillatory state, as proposed in Fischer & Spohn (1990).