Tidally Heated Exomoons around Gas Giants

Thousands of exoplanets have been discovered; however, the detection of exomoons remains elusive. Tidally heated exomoons have been proposed as candidate targets for observation; vigorous tidal dissipation can raise the moon’s surface temperature, making direct imaging possible, and cause widespread volcanism that can have a signature in transits. We assess whether the required amounts of tidal dissipation can be attained and how long it can be sustained. In a first step, we look at the thermal state of a super-Io for different orbital configurations. We show that close-in exomoons with moderate (e∼eIo ) to high ( e∼0.1 ) orbital eccentricities can feature surface heat fluxes 1–3 orders of magnitude higher than that of Io if heat transfer is dominated by heat piping or the moon has a magma ocean. In a second step, we investigate the longevity of a super-Io. The free eccentricity of an isolated close-in exomoon is quickly dampened due to tides; high orbital eccentricities can be maintained if the moon is in a mean-motion resonance with another moon and the planet is highly dissipative. However, this scenario leads to fast orbital migration. For a Mars-sized exomoon, we find that tides alone can raise the surface temperatures to more than 400 K for 10 million yr, and surface heat fluxes higher than that of Io can be maintained for billions of years. Such tidally active bodies are expected to feature more vigorous volcanic activity than Io. The material outgassed via volcanism might be detected in transits.


Introduction
With the list of confirmed exoplanets continually growing, it is no surprise that the search for exomoons has started to gain momentum. The small size of the satellites with respect to their host planets makes their detection challenging. However, with improvements in observation capabilities and novel detection methods, the observation of exomoons is now within reach. Several techniques have been proposed for the detection of exomoons (e.g., Heller 2017). These methods include transits (Kipping 2009;Ben-Jaffel & Ballester 2014;Hippke 2015;Kipping et al. 2015;Teachey et al. 2017;Teachey & Kipping 2018), microlensing (Han & Han 2002;Bennett et al. 2014), cyclotron radio emissions (Noyola et al. 2014(Noyola et al. , 2016, pulsar timing variations (Lewis et al. 2008), and direct imaging (Cabrera & Schneider 2007;Peters & Turner 2013;Heller 2016).
Tidally heated exomoons are promising targets in the exomoon hunt (e.g., Peters & Turner 2013;Ben-Jaffel & Ballester 2014;Heller 2017;Oza et al. 2019). Tidal dissipation within an exomoon can heat its interior and result in vigorous observable geologic activity. Io, the innermost Galilean satellite, is a good archetype. Although it has a radius four times smaller than that of Earth, tides raised by Jupiter result in an intrinsic surface heat flux roughly 30 times higher than that of Earth (Lainey et al. 2009;Turcotte & Schubert 2014). Tidal dissipation drives widespread volcanism, which in turn is responsible for the formation of a secondary atmosphere that extends over more than 400 Jupiter radii (Mendillo et al. 1990).
Among the outgassed material, there is sodium and potassium, which have not been detected in the gaseous envelopes of gas giants and thus can be used as a proxy for volcanic activity within a planetary system (Johnson & Huggins 2006;Oza et al. 2019). Moreover, the interaction of the outgassed secondary atmosphere with Jupiter's magnetosphere produces Io's plasma torus. While Io is unique in the solar system, objects with similar or higher levels of internal heating, super-Ios, might be common. In fact, it is possible that these kinds of objects have already been detected. Ben-Jaffel & Ballester (2014) suggested that the early ingress of close-in exoplanets WASP-12b and HD 189733b observed in the UV can be explained by the presence of a plasma torus, and Oza et al. (2019) proposed that Na signatures in the spectra of the hot Jupiter WASP-49b are evidence of a tidally heated exomoon. If confirmed, these exomoons are likely very different from Io, though; due to the close proximity of the planet to the star, surface temperatures are between 1000and3000 K, and tidal dissipation is the result of stellar tides. Transits of temperate/cold gas giants are less likely than transits of close-in giant planets, as they orbit further away from the star and thus have smaller transit probabilities (Dalba et al. 2015). However, cold gas giants suited for transit spectroscopy have already been identified, for example, HIP 41378 f (Becker et al. 2018).
In extreme cases, tidal heating can become the dominant heat source exceeding solar irradiation and have an observable footprint in the surface temperature. Peters & Turner (2013) proposed that a super-Io orbiting a cold gas giant could be directly imaged using current and planned telescopes. For this to happen, the moon should have a highly eccentric and/or short-period orbit around a cold planet. Peters & Turner (2013) concluded that Earth-sized exomoons with surface temperatures over 600 and 300 K could be detected with the Spitzer Space Telescope and future telescopes such as the James Webb Space Telescope or SPICA, respectively. The high surface temperatures needed to observe tidally heated exomoons require surface heat fluxes on the order of 500 W m −2 , 2 orders of magnitude higher than current values of Io. It remains uncertain if such extreme cases of super-Ios can exist. Furthermore, to be observable, they would have to persist for a long time, while close-in exomoons rapidly migrate, and their orbits are quickly circularized. Motivated by the prospects of detecting a super-Io around a cold gas giant, we ask the question of under which conditions can a super-Io persist. The question can be further split into (1) what is the possible thermal state of a super-Io, and (2) how long can a super-Io persist in an observable state?
In Section 2, we tackle the first question. We investigate thermal equilibrium states of exomoons of different sizes. While bigger exomoons are more prone to becoming super-Ios, current formation models limits the permissible size of exomoons around gas giants to approximately the size of Mars (Heller & Pudritz 2015a, 2015b. We push this limit and consider exomoons ranging from Io to Earth sizes and identify thermal states for which the amount of generated internal heat equals the amount of heat removed from the interior. We consider rocky exomoons with Io-like structure and composition. To account for the high temperatures reached within the mantle of a super-Io, our model allows for the formation of a sublayer of melt (Section 2.1). We compute tidal dissipation using the viscoelastic theory for self-gravitating bodies (e.g., Peltier 1974;Wu & Peltier 1982) using Andrade rheology (Section 2.2) and compare it with estimates of heat removed from the mantle (Section 2.3) via convection or melt advection (heat piping) to find thermal equilibria states. We first apply the model to Io and find that it can successfully explain its thermal state. Equipped with this model, we compute surface temperatures and heat fluxes for a range of fixed eccentricities and orbital periods to assess whether these bodies could be directly imaged or could exhibit vigorous tidal activity and have substantial exospheres.
In Section 3, we tackle the second question: how long can a super-Io live? Instead of considering fixed values for orbital period and eccentricity, we take into account the feedback between tidal dissipation, internal structure, and orbital parameters. We consider two scenarios: an isolated moonplanet system and a system with two exomoons where orbital resonances can occur. Orbital resonances are common in the solar system (e.g., Peale 1976) and responsible for the high geological activity featured on some of the outer planet moons. Classic examples of resonances in the solar system include Janus-Epimetheus, Mimas-Tethys, and Enceladus-Dione in the Saturnian system and the Laplace resonance involving the three inner Galilean moons. We consider the 2:1 mean-motion resonance (MMR) and study the thermal-orbital evolution of a Mars-sized exomoon in more detail.

Interior Structure and Rheology
We consider three different moons with radii equal to ( ) · R 1, 2, 4 Io , which represent Io and roughly Mars-and Earthsized 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 r c and radius R c and a silicate rocky layer with density r 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 sulfurrich 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). 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): 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 −1 . 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 (f) of the layer is computed as with T being the average temperature of the layer. T s and T l are the averages of the pressure-dependent solidus and liquidus temperatures (Takahashi 1990), with the pressure in gigapascals. When a slope of = dT dP 10 s l , K GPa −1 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, where h 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 −1 close to the surface to 540 KJ mol −1 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 (f) 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):

10
. 7 T T s 1 1600 2 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 m 0 . m 1 is an empirical constant, and m 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 f > 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.

Internal Heat
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 References.
(1) Anderson et al. (1996) 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) 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 ( ) I k 2 is the imaginary component of the k 2 Love number, which is related to the quality factor (Q) as . 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 0.1, and orders-of-magnitude differences can arise for very high eccentricities ( > e 0.6; Renaud et al. 2021). Additionally, for > e 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 ( ) 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˜( Y is a vector containing the Fourier-transformed radial and tangential displacements (ỹ 1 ,ỹ 2 ), the radial and shear stress (ỹ 3 , y 4 ), the gravitational potential (ỹ 5 ), and the potential stress (ỹ 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 (h m). 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 functionJ :˜( 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 z » 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 ( ) I k 2 depends on the core size, which can differ for different exomoons. However, its effect is small. For an Iosized 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.

Heat Transport
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 2001Moore , 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.

Mantle Convection
Due to the strong dependence of viscosity on temperature, most terrestrial planets in the solar system are in the stagnantlid 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 where k is the thermal conductivity of the mantle, DT 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 θ: For Newtonian viscosity, scaling arguments lead to (Solomatov 1995;Reese et al. 1999) with » a 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): For typical mantle conditions, g c is of the order of »0.01 (Nimmo & Stevenson 2000). The thickness of the elastic lithosphere (D l ) is computed as 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.

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 ( r D ), the surface gravity (g), the melt's viscosity (h l ), and a constant (τ) closely linked to the permeability exponent: 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 highefficiency melt transport scenario with n=3 and g = g g 10 6 Io and n=2 and g =g g 10 5 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.

Convection in a Magma Ocean
When the disaggregation point is reached, the sublayer of melt behaves as a magma ocean. The viscosity of a magma ocean is very low (∼0.1 Pa s), rendering heat transport via convection very efficient (Solomatov 2007). As for the tidal response, we do not model heat transport in this regime. We assume that heat transport in a magma ocean is sufficiently efficient to remove all of the internal heat.

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: Q s is the surface heat flow, and  Q 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: . 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. corresponds to a body with Andrade rheology. Stable and unstable equilibrium points are indicated in red and blue, respectively. Tidal dissipation is displayed for Io's present eccentricity and orbital frequency; a decrease in eccentricity or orbital frequency changes the amount of tidal dissipation (gray line) and drives the system evolution, as indicated by the blue line.
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 T 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.

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, »1200 K. As eccentricity increases, the equilibrium state moves to higher mantle temperatures. For moderate eccentricities (~e 2 Io ), the dislocation creep mechanism can lead to situations where   » dQ dT dQ dT s i n t . In such cases, a quasiequilibrium 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).
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 (»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 −1 . This rate is higher than the minimum 0.02 cm yr −1 required by the lack of impact craters and consistent with the 0.4-14 cm yr −1 resurfacing rate estimated from surface changes observed during Galileo's mission (Phillips 2000).

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 , 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 (  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,  M o,Io , and a tidal efficiency, h T , (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 R 2 Io and R 4 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 (  µ E 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., R 2 Io and R 4 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 R 2 Io and R 4 Io moon, respectively. Additionally, for = R R 4 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, 2 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 2 R or at its current orbital distance with an eccentricity of »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.
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?

Longevity of a Super-Io
As we have shown in the previous section, for tidally active exomoons to be observable, vigorous tidal dissipation linked to a high orbital eccentricity and/or low orbital period is needed. However, it remains to be seen if these orbital states are attainable, let alone long-lived, so that they provide opportunities for observation. Instead of taking orbital period and eccentricity as constant parameters, as we did in the first part of the paper, in this section, we explore the coupling between the interior thermal state of a moon and its orbital parameters. Clearly, we cannot tackle all possible orbital evolution scenarios; we pick two representative cases. We start by considering tidal interactions between a gas giant and an orbiting exomoon in an isolated moon-planet system. Afterward, we consider how orbital resonances can excite the orbital eccentricity of an exomoon; we focus on the simplest orbital resonance: a 2:1 MMR.

Isolated Moon-Planet System
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 (  q) are given by (Appendix C)   occurs depends on the internal structure of the moon via the second-degree Love number ( ) I k m 2, . 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  q » d dt 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, ( ) I k p 2, . In classical tidal theory, it is considered to be independent of the forcing frequency; we make this assumption and consider a range of frequencyindependent ( ) I k p 2, . However, new observations indicate a strong dependence of ( ) I k p 2, on forcing frequency (Lainey et al. 2017(Lainey et al. , 2020, which has been linked with the excitation of internal waves in the planet's gaseous envelope (e.g., Ogilvie & Lin 2004). The ( ) I k p 2, 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 ( ) I k m 2, similar to that of Io (Lainey et al. 2009). t e strongly depends on orbital distance ( 13 2 ); if the moon is far from the planet ( ⪆ a a 20 R ), the circularization timescale is more than 1 Gyr. At such orbital distances, tidal dissipation is quite low. However, a highly eccentric ( > e 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.
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.

2:1 Mean-motion Resonance
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 1 , n 2 ) are related via , where w 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 q 1, 1). In a 2:1 orbital resonance, the resonant variable =v n n 2 1 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):  Here C 1 and C 2 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 e , 1 2 is a positive number, which is small provided the eccentricities are also small. Different orbital migration timescales (t t ¹ n n 1 2 ) 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 Figure 6. t e as a function of orbital distance and moon-to-planet mass ratio for a body with | ( )| I =k 10 m 2, 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, t n and t 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)  where h/a is the aspect ratio of the disk, and M d is the disk mass. 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 (t t n n 1 2 ), 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 t 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 t n with orbital distance (Equation (24b)), this requires that | ( ( ))| | ( ( ))| I I > k n k n p p 2, 1 2, 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 , » e 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).
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 (  = v 0), and an equilibrium eccentricity is reached. The value of the equilibrium eccentricity can be computed as (Equation (29) The value of the equilibrium eccentricity is proportional to t t e 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 M 2 2 1 . The value of the equilibrium eccentricity is e 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 t e , or the system relaxes to the equilibrium configuration of the tidally driven scenario.
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, t t < 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-ofmagnitude 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 ( ) I k p 2, using Equation (29) (Figure 9(b)). Moderate-to-high orbital eccentricities ( ---10 10 3 1 ) that result in high surface heat fluxes are attained, provided that (1) the moon-to-planet mass ratio is high 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 thermalorbital 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 a 2 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 ( ( ) I k m 2, ), 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 | ( )| I k p 2, is frequency-independent, and we vary its value from 10 −3 to 10 −6 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
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 | ( )| I k p 2, , with higher values leading to higher values of orbital eccentricity. Meanwhile, the moon migrates at a rate dependent on ( ) I k p 2, , tidal dissipation decreases, and the moon cools down enough to enter the heat-pipe regime. Upon cooling, the value of ( ) I k m 2, changes, and, in a range of 1-10 Myr, a new equilibrium orbital eccentricity is attained.
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 | ( )| I =k 10 p 2, 4 , surface heat fluxes of more than 1000 W m −2 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 heatpipe 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 (t 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 | ( )| I = -k 10 p 2, 5, 6 is lower than that for | ( )| I = -k 10 p 2, 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).

Conclusions
We started this paper with a clear question-what are the prospects of detecting a super-Io?-and we addressed two important sides of it: where and in which thermal state we can expect to find super-Ios and for how long we expect a super-Io to be tidally active.
To do so, we presented a thermal model of an exomoon. Based on our current knowledge of Io, we considered a multilayered model that allows for the formation of a sublayer of melt in which heat can be transported via melt advection. Our model confirms the findings of Moore (2001Moore ( , 2003 and Moore & Webb (2013) that this mechanism plays a crucial role in the heat budget of a tidally heated (exo)moon. We applied our model to Io and found that it can successfully explain its thermal state and is consistent with Galileo's observation hinting at a partially melted asthenosphere (Khurana et al. 2011).
In order to simplify our model, we made some assumptions. As is commonly done (e.g., Fischer & Spohn 1990;Hussmann & Spohn 2004;Henning et al. 2009;Shoji & Kurita 2014;Renaud & Henning 2018), we used the viscoelastic tidal theory to model tidal dissipation. However, when tidal dissipation is very high, the melt fraction increases and can reach high values. The body is better described either as a highly porous material consisting of a matrix of rock filled with magma or, for higher melt fractions, as a magma ocean with rocky particles in suspension. Our treatment of tidal dissipation and heat transport in the high-melt fraction regime was highly simplified. In such circumstances, the classical theory of tides in solid viscoelastic does not adequately describe the behavior of the body. Advances have been made in understanding the role of porosity (Liao et al. 2020) and tides in a liquid reservoir (e.g., Tyler et al. 2015;Matsuyama et al. 2018;Rovira-Navarro et al. 2019;Hay et al. 2020) in tidal dissipation, yet our knowledge remains limited. Further study of these regimes will be key to understanding Io and ultimately super-Ios. Additionally, we considered that material properties only vary radially and not laterally. While this assumption allows a first-order understanding of Io and super-Ios, unveiling their complex dynamics might require taking into account the complex feedback between tidal dissipation, internal properties, and resulting heterogeneities (Steinke et al. 2020). In the case of Io, the availability of high-resolution observations justifies this kind of modeling, but we are still very far from this point for exomoons. Figure 10. Thermal-orbital evolution of a Mars-sized exomoon orbiting around a 5 M ♃ planet in a 2:1 orbital resonance with an exomoon of the same size. Panel (a) shows the semimajor axis in terms of the Roche limit (y-axis) as a function of time (x-axis). The trajectories are colored according to surface heat flux and temperature. Panel (b) shows the detail of the evolution of the eccentricity (solid lines) and semimajor axis (dashed lines), and panel (c) shows the surface temperature using a logarithm timescale. We applied our model to exomoons ranging from Io to Earth-sized and found that exomoons orbiting within a few Roche radii of their planet experience high levels of tidal dissipation if they have moderate eccentricities (similar to Io). Thermal equilibrium states with surface heat fluxes higher than those of Io are obtained for a wide range of orbital parameters. For some of these cases, tidal dissipation is very high, and the thermal budget of the moon is dominated by tidal dissipation, resulting in a high surface temperature that could make direct imaging by future telescopes possible (Peters & Turner 2013).
While our analysis of thermal equilibrium states as a function of orbital parameters showed that super-Ios can result for different orbital configurations, one key question remained: do we expect a population of exomoons to fill this space, and for how long? This question should be tackled from a moon formation and evolution point of view. We addressed the second aspect by considering the coupling between interior structure and orbit evolution of a close-in exomoon.
We showed that an Earth-sized exomoon orbiting far away from the planet ( / > a a 20 R ) with a high eccentricity can remain tidally active for billions of years and speculated that this might happen if such an exomoon is captured via a mechanism known as the binary-exchange capture process (Williams 2013). However, under these circumstances, our assumption of the moon being tidally locked and having a small obliquity could break. Moreover, higher-order eccentricity terms will be become relevant for the orbital evolution of the moon (Renaud et al. 2021). The more general problem of capture and tides for a nonsynchronous rotating moon should be addressed. This opens the door to higherorder spin-orbit resonances, such as a Mercury-type 3: 2 resonance (e.g., Dobrovolskis 2007;Makarov 2012;Walterová & Běhounková 2020), as well as obliquity tides that can be investigated in future work.
For closer exomoons, circularization timescales are on the order of million of years, which implies that long-lived super-Ios require active forcing of the eccentricity, for example, via MMRs (e.g., Peale 1976). We considered a super-Io in a 2: 1 orbital resonance and found a trade-off between high orbital eccentricity and fast orbital migration. We obtained that in order for the forced eccentricity to be high enough for surface temperatures to be dominated by tidal heating, the orbited exoplanet should be highly dissipative (i.e., have a low quality factor). Low quality factors are not compatible with the equilibrium tide theory of Goldreich & Nicholson (1977) but are in agreement with the low values measured for Jupiter and Saturn (Lainey et al. 2009(Lainey et al. , 2020. However, we showed that highly dissipative planets lead to fast orbital migration, which limits the amount of time the moon spends near the planet, where tidal dissipation is high. We studied in more detail the case of a Mars-sized exomoon orbiting a super-Jovian planet. We found that tidal dissipation decreases below solar irradiation »500 million yr after the establishing of the MMR, limiting the window for direct imaging. However, the moon can remain tidally active for billion of years with heat fluxes exceeding those of Io. This implies that a plasma torus denser than that of Jupiter can be expected in the system, which provides an opportunity for detection in transmission spectra (Johnson & Huggins 2006;Oza et al. 2019). Finding a tidally heated exomoon via transit spectroscopy or direct imaging would provide an important constraint on the planet's interior structure due to the relation between the moon's eccentricity and dissipation within the planet. The solution to the previous set of differential equations can be written as˜( where ɷ is the forcing frequency, and Y is the so-called fundamental matrix, which equals and C i is a vector of integration constants. The solution at the surface of the moon (˜( ) y R ) can be computed by propagating the solution from the core-mantle boundary (R c ) to the surface (R) by imposing continuity at the layers' boundaries: 1 1 I c is a matrix that follows from the boundary conditions at the core-mantle boundary and is given by and C c is a vector of three integration constants that we obtain by applying the surface boundary conditions: For different values of volumetric heat (q vol ), we integrate the previous set of equations from the bottom of the asthenosphere to the top. The average melt fraction f is then computed, and a curve relating the average melt fraction (f) and q vol is obtained. The effects of planet and moon tides on the evolution of a moon-planet system have been widely studied (e.g., Kaula 1964;Boué & Efroimsky 2019). Here we follow Efroimsky & Williams (2009) and Boué & Efroimsky (2019) to obtain the tidal effects on the orbit of a synchronously rotating satellite with low inclination orbiting around a planet spinning at angular rate  q. The evolution of the satellite's orbital frequency and eccentricity are given by Boué & Efroimsky (2019, Equations (143) and (156)), and the change of the planet's spin rate due to the moon is given by Efroimsky & Williams (2009, Equation (34)): Finally, it is important to note that higher-order terms can become important for high orbital eccentricities ( > e 0.1; Renaud et al. 2021).
Comparing Equation (C4) with the equations used in previous studies (e.g. , Yoder 1979;Malhotra 1991;Fischer & Spohn 1990;Shoji & Kurita 2014), we note that those studies use p e = 3 instead of the value obtained above. This discrepancy arises from the incorrect assumption in some past publications that eccentricity damping occurs at constant angular momentum. The thermal-orbital evolution scenarios presented in Section 3 are markedly different from those obtained by Fischer & Spohn (1990) and Hussmann & Spohn (2004) for Io. In Fischer & Spohn (1990) and Hussmann & Spohn (2004), the thermal-orbital evolution is characterized by a nearly equilibrium phase followed by an oscillatory phase of alternating cold and warm phases and finally a runaway cooling of the body. In contrast, in our model, the body evolves following a series of near-equilibrium thermal-orbital states. The difference is due to the use of Andrade rheology instead of Maxwell rheology, as well as the introduction of a more efficient heat transport mechanism: heat piping. Here we briefly explore how these two factors affect the thermal-orbital evolution of a rocky exomoon.
The difference can be explained in terms of the location and stability of thermal equilibrium points. To illustrate this point, we consider the thermal-orbital evolution of an Io-sized moon using the two different models presented in Section 2.4: (a) Maxwell rheology and heat transport via mantle convection and (b) Andrade rheology and heat transport via heat piping. In both cases, we consider that the moon starts its evolution at = a a 2 R with an initial orbital eccentricity of 10 −3 in an orbital resonance with a moon of the same size. Due to the close proximity of the moon to the planet, tidal dissipation is high, and a stable equilibrium point with > T T res is reached in both cases. As the moon migrates outward and tidal dissipation decreases, the equilibrium mantle temperature decreases accordingly until the point = T T res . This point is unstable for model (a) but stable for model (b). For model (a), further orbital migration starts a runaway cooling phase. As the moon cools downs, ( ) I k m 2 sharply decreases, which leads to an increase of the orbital eccentricity (Equation (29)) and thus tidal dissipation, causing the body to heat up again (Figure 11). This process is repeated several times, resulting in thermal-orbital oscillations. In contrast, when Andrade rheology and heat piping are included, the equilibrium point is stable and the moon evolves following thermal equilibrium states, no oscillatory phase occurs.
This result can also be interpreted in the context of the linear stability analysis presented in Ojakangas & Stevenson (1986). Ojakangas & Stevenson (1986) showed that the stability and subsequent orbit evolution depends on the exponents of the power dependence of tidal dissipation and heat transport with temperature (  = n d Q dT ln int and  = m d Q dT ln s ). The Andrade rheology and heat piping reduce n and increase m, which brings the system to the stable regime. Figure 11. Thermal-orbital evolution of an Io-sized exomoon in a 2:1 orbital resonance around a Jovian planet. The solid line corresponds to a model with Maxwell rheology and heat transport dominated by convection, and the dashed line corresponds to Andrade rheology and heat transport dominated by heat piping. In both cases, | ( ) | I =k 10 p 2 5 is assumed.