# The tidal-thermal evolution of the Pluto-Charon system

Amirhossein Bagherian, Amir Khan, Frédéric Deschamps, Henri Samuel, Mikhail Kruglyakov, Domenico Giardini

► To cite this version:

Amirhossein Bagherian, Amir Khan, Frédéric Deschamps, Henri Samuel, Mikhail Kruglyakov, et al.. The tidal-thermal evolution of the Pluto-Charon system. Icarus, 2022, 376, 10.1016/j.icarus.2021.114871. insu-03643024

HAL Id: insu-03643024

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

Submitted on 20 Apr 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research Paper

# The tidal-thermal evolution of the Pluto-Charon system

Amirhossein Bagheri<sup>a,*</sup>, Amir Khan<sup>a,b</sup>, Frédéric Deschamps<sup>c</sup>, Henri Samuel<sup>d</sup>,
Mikhail Kruglyakov<sup>e,f</sup>, Domenico Giardini<sup>a</sup>

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

## ARTICLE INFO

**Keywords:**

Pluto

Charon

Trans-Neptunian objects

Tidal evolution

Thermal evolution

Subsurface ocean

## ABSTRACT

The existence of subsurface oceans on the satellites of the giant planets and Trans-Neptunian objects has been predicted for some time. Liquid oceans on icy worlds, if present, exert a considerable influence on the dynamics of the ice-ocean system and, because of the astrobiological potential, represent an important objective for future missions to the outer solar system. The Pluto-Charon system is representative of an icy moon orbiting a dwarf planet that is believed to have formed from the remnants of a giant impact. The evolution of icy moons is primarily controlled by the mode and efficiency of heat transfer through the outer ice shell, which is influenced by the presence of impurities, by tidal dissipation in the ice shell, and by the radioactive element budget in the silicate core. Previous studies on the evolution of the Pluto-Charon system generally considered either only the thermal or the tidal evolution, and in the cases where both were considered, the important effect of the presence of impurities in the liquid oceans was not addressed. Here, we consider the joint tidal-thermal evolution of the Pluto-Charon system by combining a comprehensive tidal model that incorporates a viscoelastic description of the tidal response with a parameterized thermal convection model developed for icy worlds. This approach enables an extensive analysis of the conditions required for the formation and maintenance of subsurface liquid oceans up to the present. Our results show that because of relatively fast circularization and synchronization of the orbits of Pluto and Charon, tidal heating is only important during the early stages of evolution (<1 Myr). As part of our study, we test the sensitivity of our results to a number of parameters that pertain to the orbital and thermal history. In all the studied cases, oceans on Pluto are always predicted to remain liquid to the present, ranging in thickness from 40 km to 150-km, whereas oceans on Charon, while in-place for approximately 4 Gyr, have solidified. This is supported by New Horizons observations of primarily extensional faults on Pluto and both extensional and compressional faults on Charon.

## 1. Introduction

The Kuiper belt harbors numerous planetary objects of diverse internal structure and surface features, including the dwarf planet Pluto. Pluto hosts five known satellites named Charon, Kerberos, Hydra, Nix, and Styx, of which Charon is by far the largest. The Pluto system is, because of Pluto's size and relative brightness, presently the best-studied of all of the Trans-Neptunian objects (TNOs) (Hussmann et al., 2006). This is a consequence of a protracted history of Earth-based remote sensing (Malhotra and Williams, 1997; Dobrovolskis et al., 1997; Olkin et al., 2003) and not least the flyby of the New Horizons spacecraft in 2015 (Spencer et al., 2020). Pluto was found to display both a complex and an active geology that encompasses an extensive

range of surface ages and a dynamic linkage between surface and atmosphere (Nimmo et al., 2017; Spencer et al., 2020). In comparison, Charon appears to be both geologically and compositionally distinct to Pluto (Spencer et al., 2021). The Pluto-Charon mass ratio was also observed to be much higher than that of the Earth-Moon system, which is untypical of Planet-satellite systems in our Solar System. Charon, in analogy with the Earth-Moon system, is believed to have formed as a result of a collision between Pluto and a Kuiper belt object (Canup, 2010; Sekine et al., 2017; Arakawa et al., 2019). Like Earth's Moon, Charon would initially have been closer to Pluto, but because of tidal dissipation within the two bodies, Charon would have been driven

<footnote><sup>a</sup> Institute of Geophysics, ETH Zürich, Zürich, Switzerland</footnote>

<footnote><sup>b</sup> Physik-Institut, University of Zürich, Zürich, Switzerland</footnote>

<footnote><sup>c</sup> Institute of Earth Sciences, Academia Sinica, 128 Academia Road, Sector 2, Nangang, Taipei 11529, Taiwan</footnote>

<footnote><sup>d</sup> Institut de Physique du Globe de Paris, CNRS, Université de Paris, Paris, France</footnote>

<footnote><sup>e</sup> University of Otago, Dunedin, New Zealand</footnote>

<footnote><sup>f</sup> Geoelectromagnetic Research Center, Institute of Physics of the Earth, Moscow, Russia</footnote>

<footnote>* Corresponding author.
E-mail address: amirhossein.bagheri@erdw.ethz.ch (A. Bagheri).</footnote>
further away until it reached its current synchronous state (Dobrovol-
skis et al., 1997).

Generally speaking, tidal evolution drives planetary systems toward
equilibrium states by damping their orbital eccentricity and forcing
the spin rates toward stable synchronous rotation, while adjusting the
separation between the planetary objects. During the deceleration of
the spin of the planetary objects, heat is produced by friction, which,
in addition to radiogenic heating, changes the thermal structure of
the planet. As the planetary objects thermally evolve, their interior
properties change, which in turn, impacts their tidal response (e.g.,
Robuchon and Nimmo, 2011; Saxena et al., 2018; Samuel et al., 2019;
Bagheri et al., 2021; Renaud et al., 2021). Hence, the tidal and ther-
mal evolution co-modulate, necessitating their joint consideration in
evaluating the evolution of planetary systems. Because the amount of
tidal dissipation depends on the thermal state, physical structure, and
orbital parameters (distance, eccentricity, and spin and orbit rates) of
the objects, these all need to be studied within a single framework.

The New Horizons mission found indications that Pluto harbors
and Charon had a subsurface ocean beneath their ice-covered sur-
faces (Nimmo et al., 2016, 2017; Olkin et al., 2017). The possibility
that both bodies are able to sustain a liquid ocean beneath an icy cover
enhances their status as prime targets for the search for extra-terrestrial
life (Vance et al., 2018). For an ocean to form and remain liquid, the
presence of long-lived heat sources is required, of which radioactive
and tidal heating are the most significant contributors (McKinnon et al.,
1997; Hussmann and Spohn, 2004; Schubert et al., 2010; Robuchon and
Nimmo, 2011; Saxena et al., 2018). Consequently, understanding the
long-term thermal and tidal evolution of planetary systems is a central
tenet in evaluating the possibility for the existence of a present-day
subsurface ocean and, in turn, the astrobiological potential of the outer
Solar System (McKinnon, 2006; Mottl et al., 2007; Vance et al., 2007).

Several studies have addressed the evolution of the Pluto–Charon
system in the context of coupled thermal-orbital evolution models (e.g.,
Robuchon and Nimmo, 2011; Barr and Collins, 2015; Hammond et al.,
2016; Desch and Neveu, 2017; Saxena et al., 2018). The role of tidal
heating in the evolution of Kuiper belt Objects, including Pluto–Charon,
was found to be comparable to and even higher than the heat produced
by the radioactive decay of long-lived isotopes. Saxena et al. (2018), for
example, observed that subsurface oceans containing a small amount of
impurities and tidal heating due to initially high spin rates may enable
liquid water and cryovolcanism to persist until the present (Moore
et al., 2016; Neveu et al., 2015; Beyer et al., 2019). Yet, these studies
generally did not consider dissipation in both bodies and relied on tidal
evolution models that are inadequate for the case of a highly eccen-
tric and/or non-synchronously rotating system (Bagheri et al., 2021;
Renaud et al., 2021; Arakawa et al., 2021). Specifically, tidal models
that truncate eccentricity functions to e² (where e is eccentricity) on
eccentric orbits (>0.1), can result in erroneous spin rate evolution,
spin–orbit resonances, and heating rates in comparison to computa-
tions that include higher-order eccentricity terms. Moreover, while
the aforementioned studies found the evolution to be fast (≤1 Myr),
incorporation of non-synchronous rotation can slow the evolution down
if higher-order spin–orbit resonances are encountered (e.g., Saxena
et al., 2018).

With this in mind, it is the purpose here to build upon and extend
earlier studies on the evolution of the Pluto–Charon system by combin-
ing a comprehensive tidal model that incorporates a proper viscoelastic
description of the tidal response of planetary bodies (Bagheri et al.,
2021) with a thermal evolution model that tracks the temperature
change in the interior over ~4.5 Gyr. Our parameterized thermal
convection model for the icy crust is based on the work of Deschamps
and Vilella (2021) and takes into account the effect of impurities in
the liquid ocean, while for the silicate core we rely on the thermal
evolution model of Samuel et al. (2019). This approach will enable
an extensive analysis of (1) the tidal–thermal evolution of the Pluto–
Charon system, allowing us to assess the relative role of radiogenic to

tidal heating and (2) the conditions required for the formation and
maintenance of subsurface liquid oceans up until the present. More
generally, the methodology and results presented here can be exploited
to understand the evolution of icy satellites and easily be extended to
TNOs and exoplanets.

The manuscript is arranged as follows. In Section 2, we describe
current observations and the constraints they provide on the interior
structure of Pluto and Charon; in Section 3, we detail the thermal
evolution and tidal models; in Section 4, we present and discuss the
results with particular emphasis on model parameter exploration and
its effect on the outcome; and summarize and conclude in Section 5.

2. Pluto and charon

2.1. Interior

Direct information bearing on the interior structure of Pluto and
Charon is, as noted, scarce. Based on New Horizons observations, mass
and radius of the two objects could be determined, allowing for the
estimation of their mean densities (Nimmo et al., 2017) (see Table 1).
The bulk densities of Pluto and Charon are clearly higher than that of
ice and therefore the interiors of the two bodies must be composed of
denser material such as silicates and possibly also metals (McKinnon
et al., 2008). The mean density of Pluto indicates a rock fraction of
about 2/3, the rest possibly consisting of ice (McKinnon et al., 2017),
whereas the lower density of Charon suggests a slightly lower rock mass
fraction. While porosity affects density, it is not expected to be able to
account for the density difference between the two bodies (McKinnon
et al., 2017; Bierson et al., 2018).

Pluto is believed to have formed from the hydrated silicate cores and icy material of the mantles of two impacting objects, both of which were already differentiated (Desch, 2015) or partially differ- entiated (Canup, 2010; Desch and Neveu, 2017). Despite the lack of measurements of the moments of inertia, New Horizons observations indicate that both Pluto and Charon are most probably differenti- ated (Stern et al., 2015, 2018; Spencer et al., 2021). The observations include surface spectra that are dominated by ices and the lack of compression in the surface geological record. Were the interiors un- differentiated, (1) surface spectra would have been expected to be contaminated by the presence of silicates and (2) the formation of deeper and denser high-pressure ices should have resulted in compres- sion as Pluto and Charon cooled (McKinnon et al., 2017; Grundy et al., 2016; Hammond et al., 2016). Instead, the observations indicate the presence of extensional tectonic features associated with the expan- sion that occurred when an early ocean froze above a differentiated interior (Stern et al., 2015; Moore et al., 2016; Beyer et al., 2017).

Based on previous estimates, radioactive heat production within
Pluto would be sufficient to melt a conductive ice shell and maintain a
global subsurface ocean until the present day (Hammond et al., 2016;
Bierson et al., 2018). In contrast, a convective ice shell would allow for
rapid removal of heat such that an ocean would never develop (Robu-
chon and Nimmo, 2011). Thus, the presence of an ocean today would
suggest a cold, rigid, and conductive ice shell that allows for little
interaction between ocean and surface. By the same argument, Charon
is not expected to host a subsurface ocean today, although the latter
might have developed at some point earlier in its evolution (Desch
and Neveu, 2017; Bierson et al., 2018). Imaged extensional tectonic
features are evidence of strains (Beyer et al., 2017) that are expected if
an ocean refreezes (Spencer et al., 2020). If Charon's initial orbit was
non-circular, it would also have experienced an early episode of tidal
heating and stress generation (Rhoden et al., 2015).

Plausible present-day interior structure models for Pluto and Charon
are shown as cross-sections in Fig. 1 (Stern et al., 2015; Bierson et al.,
2018; McKinnon et al., 2017; Nimmo et al., 2017; Rhoden et al.,
2020). While a small iron core is indicated in the cross-section, we
consider the core to belong to the silicate part hereafter, because of
<table><caption>**Table 1**
Properties of the Pluto–Charon system are taken from Nimmo et al. (2017) Stern et al. (2015), and Brozović et al. (2015).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Pluto</th><th>Charon</th></tr></thead><tbody><tr><td>Radius (m)</td><td><i>R</i></td><td>1.1883 × 10<sup>6</sup></td><td>0.606 × 10<sup>6</sup></td></tr><tr><td>Mass (kg)</td><td><i>M</i></td><td>1.328 × 10<sup>22</sup></td><td>1.603 × 10<sup>21</sup></td></tr><tr><td>Density (g/cm<sup>3</sup>)</td><td><i>ρ</i></td><td>1.854</td><td>1.701</td></tr><tr><td>Semi-major axis (m)</td><td><i>a</i></td><td>—</td><td>19.596 × 10<sup>6</sup></td></tr><tr><td>Spin period (day)</td><td><i>θ</i></td><td>6.387</td><td>6.387</td></tr><tr><td>Eccentricity</td><td><i>e</i></td><td>—</td><td>0.00005</td></tr><tr><td>Inclination</td><td><i>i</i></td><td>—</td><td>0.0</td></tr></tbody></table>

<table><caption>**Table 2**
Model parameter values for the reference model. Values for ice and water are computed using SeaFreeze (Journaux et al., 2020) at the surface temperature of Pluto, while those for silicate are taken from Jackson and Faul (2010) and McKinnon et al. (2017).</caption><thead><tr><th>Parameter/unit</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td><strong>Ice</strong></td><td></td><td></td></tr><tr><td>Shear modulus (GPa)</td><td>μ<sub>ice</sub></td><td>4.8</td></tr><tr><td>Bulk modulus (GPa)</td><td>κ<sub>ice</sub></td><td>10.5</td></tr><tr><td>Density (g/cm<sup>3</sup>)</td><td>ρ<sub>ice</sub></td><td>0.92</td></tr><tr><td><strong>Water</strong></td><td></td><td></td></tr><tr><td>Shear modulus (GPa)</td><td>μ<sub>w</sub></td><td>0</td></tr><tr><td>Bulk modulus (GPa)</td><td>κ<sub>w</sub></td><td>2.2</td></tr><tr><td>Density (g/cm<sup>3</sup>)</td><td>ρ<sub>w</sub></td><td>1.0</td></tr><tr><td><strong>Silicate</strong></td><td></td><td></td></tr><tr><td>Shear modulus (GPa)</td><td>μ<sub>sil</sub></td><td>60</td></tr><tr><td>Bulk modulus (GPa)</td><td>κ<sub>sil</sub></td><td>120</td></tr><tr><td>Density (g/cm<sup>3</sup>)</td><td>ρ<sub>sil</sub></td><td>3.0–3.6</td></tr></tbody></table>

the negligible effect of a tiny core on the evolution. Given the dearth of direct observations on the interior, the radial extent of each layer is not well-constrained, which is reflected in the thickness ranges indicated in the cross-sections. Estimations of the density and volume of the silicate part and ice/water layer are such that the measured mean densities for Pluto and Charon (Table 1) are satisfied. For the models of Pluto and Charon adopted here, we vary the density of the silicate core between 3–3.6 g/cm³, relying on the studies of McCord and Sotin (2005) and Feistel and Wagner (2006), as a consequence of which the radius of the silicate part varies in the range 820–890 km for Pluto and 395–430 km for Charon, to satisfy the mean density of the bodies. In turn, this implies ice/water thicknesses of 295–370 km and 175–210 km for Pluto and Charon, respectively. Physical properties for each layer are compiled in Table 2. The elastic values of the ice layer are computed at the surface temperature of Pluto using SeaFreeze (Journaux et al., 2020) and are therefore slightly higher than the average elastic modulus of the ice shell.

## 2.2. Orbital properties

The observed orbital properties, including semi-major axis (*a*), eccentricity (*e*), inclination (*i*), and spin rate (*θ*) of Pluto and Charon are compiled in Table 1. The two bodies are presently tidally locked to each other, i.e., their spin periods are equal to their orbital period around their center of mass, ensuring that each always presents the same face to the other, and their respective orbits about the center-of-mass are almost perfectly circular. The present-day eccentricity and inclination of the orbit is very close to zero (Buie et al., 2012; Stern et al., 2015; Nimmo et al., 2017; Brozović et al., 2015). This implies that there are no significant tidal forces acting on the two bodies, and hence, negligible tidal dissipation occurring within either object at the present.

# 3. Methods

## 3.1. Thermal evolution of Pluto and Charon

The thermal evolution of a planetary object is controlled by the processes that produce heat and those responsible for transferring the heat within the different planetary envelopes to the surface, either through radiation, conduction or convection. The main heat sources in a binary system are:

1. the impact heating associated with accretion during planet formation,

2. the gravitational energy released during planetary differentiation,

3. radiogenic heating in the silicate component from the decay of long-lived radioactive isotopes (U, Th, and K), and

4. tidal heating due to viscoelastic dissipation.

Of these sources, only (3) and (4) are of relevance for the long-term evolution of the planet. Sources related to (1) and (2) are principally linked to the very earliest stages of planetary accretion. A large part of the accretional heat is likely to have radiated immediately into space (Hussmann et al., 2010), whereas the energy release during early differentiation is estimated to be 10% at most of the accretional energy (Schubert et al., 1986) and consequently, negligible in comparison to the contribution from e.g., radiogenic heating (Hussmann and Spohn, 2004). Since we consider the evolution of the Pluto–Charon system from the time postdating formation and differentiation, we neglect their contributions in the following. In summary, the main heat sources affecting the thermal evolution of Pluto and Charon are those associated with radiogenic heating in the silicate parts and tidal dissipation in the ice shell. Tidal dissipation in the solid core is negligible because of the low temperatures (see Section 3.2.2). Also, since we are assuming the cores of Pluto and Charon to be solid with little porosity, any associated tidal dissipation (cf. Choblet et al., 2017) is considered to be negligible.

As illustrated in Fig. 1c, the evolution of the ice shell is computed based on the energy balance between (1) the heat produced in the silicate part (radiogenics) and entering the ice shell from below; (2) the heat produced in the shell (tides and crystallization/melting of ice); and (3) the heat that can be transported from within the ice shell to the surface. Any ocean that is either initially present or forms during the course of the evolution, is assumed to be adiabatic. We also assume that any radiogenic heating leached from the core as a result of hydrothermal activity (Castillo-Rogez and Lunine, 2010) and dissipation from ocean currents (Tyler, 2008) are negligible. A detailed treatment of these effects is beyond the scope of this study.

The properties of the outer ice shells of Pluto and Charon that are made up of ice Ih may allow thermal convection to operate within them (Hussmann et al., 2007). Details of this mechanism depend on the physical properties of the system. The rheology of ice plays an important role in heat transfer through the outer ice shell and depends strongly on the material temperature. Because the temperature-dependence of the viscosity of ice Ih is high (e.g., Durham et al., 2010), convection in the outer ice shell may occur in the so-called stagnant-lid regime, which reduces the amount of heat that can be transported to the surface because of the presence of a thermally conductive and rigid (high viscosity) lid at the top of the system (e.g., Moresi and Solomatov, 1995; Davaille and Jaupart, 1993). Stagnant-lid convection strongly influences the heat flux and interior temperature, which both depend on the top-to-bottom thermal viscosity contrast of the ice layer. The release of heat through tidal dissipation (to be discussed in Section 3.2) within the ice shell further influences the properties of the ice shell and its ability to transport heat to the surface, as suggested by numerical simulations of mixed-heated (basal and internal) thermal convection (e.g., Travis and Olson, 1994; Deschamps et al., 2010).

In the ice and water layers, we model thermal evolution using the parameterized convection model of Grasset and Sotin (1996) modified
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Cross-sections of plausible interior structure models for Charon (a) and Pluto (b). In (a) and (b), the thicknesses of the silicate parts (and therefore also of the ice layers) range from 395–430 km for Charon and 820–890 km for Pluto, respectively, to keep the mean densities (Table 1) of the bodies fixed. Attenuation values (Q) are also indicated for the various layers. (c) illustrates the location of the main heat sources that are considered in this study, while arrows indicate the direction of the heat flow.</figcaption></figure>

by (Deschamps and Vilella, 2021) to account for the effects of mixed-heated stagnant-lid convection within the outer ice shell. This approach uses scaling laws derived from simulations of thermal convection to estimate the average temperature within the ice layer and the heat flux at its bottom. The growth of the ice shell (thickness) is then estimated from the difference between the heat flux coming from the silicate part (hereinafter core) and the heat flux entering the ice shell. At the boundary between the ice shell and the subsurface ocean, energy conservation is written

$$ \frac{dr_{bot}}{dt} \left[ \rho_w C_w \left( \frac{\partial T_{bot}}{\partial r} - \frac{\partial T_{ad}}{\partial r} \right) \frac{r_{bot}^3 - r_c^3}{3} - \rho_I L_I r_{bot}^2 \right] = r_{bot}^2 \phi_{bot} - r_c^2 \phi_c, \quad (1) $$

where *t* is time, $T_{bot}$ is the temperature at the bottom of the ice layer given by the liquidus of the ocean, $\phi_{bot}$ is the heat flux at the bottom of the ice layer, $r_{bot}$ is the radius of the bottom of the ice layer, $r_c$ is the core radius, $\phi_c$ is the heat flux at the top of the core, $\rho_w$ and $C_w$ are the liquid water density and heat capacity, respectively, and $\rho_I$ and $L_I$ the density and latent heat of fusion of ice Ih, respectively. $T_{ad}$ is the adiabatic temperature in the ocean and is given by

$$ T_{ad}(r) = T_{bot}(r_{bot}) \left[ 1 - \frac{\alpha_w}{\rho_w C_w} \rho_I g (r - r_{bot}) \right], \quad (2) $$

where $\alpha_w$ is the coefficient of thermal expansion of liquid water, *g* is gravity, $\rho_I$ is the density of ice Ih, and $r_{top}$ is the radius of the top of the ocean layer. We model the ocean as an inviscid fluid layer through which heat is transported immediately from the bottom to the top of the ocean. We assume a uniform temperature profile for the ocean and ignore the tiny increase in temperature resulting from the adiabatic temperature gradient. Solving Eq. (1) for $r_{bot}$, i.e., ice layer thickness, requires knowledge of the heat fluxes at the top of the core ($\phi_c$) and at the bottom of the ice shell ($\phi_{bot}$) with time.

The heat flux at the top of the core is obtained by modeling its thermal evolution, which is governed by solving the time-dependent heat diffusion equation following the approach of Samuel et al. (2019). We assume that the core has a carbonaceous chondrite composition (Lodders, 2003) and consider four radioactive elements: $^{235}$U, $^{238}$U, $^{232}$Th, and $^{40}$K. The averaged radial temperature profile of the core ($T_c$) and the heat flux at its top ($\phi_c$) are calculated by solving the heat diffusion equation using

$$ \rho_c C_c \frac{\partial T_c}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 k_c \frac{\partial T_c}{\partial r} \right) + H(t), \quad (3) $$

where $\rho_c$, $C_c$, and $k_c$ are the density, the specific heat, and thermal conductivity of the core, respectively, and the volumetric internal heating rate, $H$, is the heat caused by tidal dissipation. The initial conditions are the initial thermal profiles for either body (shown in Fig. 3). Values for all parameters are listed in Table 3.

<caption>**Table 3**</caption>

Thermal properties of the different layers of the model (after Deschamps and Vilella (2021) for the icy crust and liquid water parts and Samuel et al. (2019) for the silicate core).

<table><thead><tr><th>Parameter (unit)</th><th>Symbol</th><th>Value/expression</th></tr></thead><tbody><tr><td><strong>Ice Ih properties</strong></td><td></td><td></td></tr><tr><td>Surface temperature K</td><td><i>T</i><sub>surf</sub></td><td>40</td></tr><tr><td>Thermal expansion (1/K)</td><td><i>α</i><sub>I</sub></td><td>1.56 × 10<sup>-4</sup></td></tr><tr><td>Thermal conductivity (W/m/K)</td><td><i>k</i></td><td>2.6</td></tr><tr><td>Heat capacity (J/kg)</td><td><i>C</i><sub>p</sub></td><td>7.037T + 185</td></tr><tr><td>Thermal diffusivity (m<sup>2</sup>/s)</td><td><i>κ</i><sub>I</sub></td><td><i>k</i>/(ρ<sub>I</sub><i>C</i><sub>p</sub>)</td></tr><tr><td>Latent heat of fusion (kJ/kg)</td><td><i>L</i><sub>I</sub></td><td>284</td></tr><tr><td>Reference bulk viscosity (Pa s)</td><td><i>η</i><sub>ref</sub></td><td>10<sup>14</sup></td></tr><tr><td>Activation energy (kJ/mol)</td><td><i>E</i><sup>*</sup></td><td>60</td></tr><tr><td>Initial thickness (Pluto) (km)</td><td><i>D</i><sup>ice</sup><sub>init</sub></td><td>298</td></tr><tr><td>Initial thickness (Charon) (km)</td><td><i>D</i><sup>ice</sup><sub>init</sub></td><td>171</td></tr><tr><td><strong>Liquid water properties</strong></td><td></td><td></td></tr><tr><td>Initial ammonia content (wt%)</td><td><i>A</i><sub>init</sub></td><td>3</td></tr><tr><td>Thermal expansion (1/K)</td><td><i>α</i><sub>w</sub></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>Heat capacity (J/kg)</td><td><i>C</i><sub>w</sub></td><td>4180</td></tr><tr><td>Initial thickness (Pluto) (km)</td><td><i>D</i><sup>w</sup><sub>init</sub></td><td>15</td></tr><tr><td>Initial thickness (Charon) (km)</td><td><i>D</i><sup>w</sup><sub>init</sub></td><td>15</td></tr><tr><td><strong>Silicate core properties</strong></td><td></td><td></td></tr><tr><td>Thermal conductivity (W/m/K)</td><td><i>k</i><sub>c</sub></td><td>2.4</td></tr><tr><td>Heat capacity (J/kg/K)</td><td><i>C</i><sub>c</sub></td><td>1100</td></tr><tr><td>Core radius (Pluto) (km)</td><td><i>R</i><sub>c</sub></td><td>875</td></tr><tr><td>Core radius (Charon) (km)</td><td><i>R</i><sub>c</sub></td><td>420</td></tr><tr><td>Density (kg/m<sup>3</sup>)</td><td><i>ρ</i><sub>c</sub></td><td>3100</td></tr></tbody></table>

Depending on its properties, principally viscosity and thickness, the ice shell may transport heat either by conduction or by convection, leading to two different estimates of heat flux at its bottom. For a static, thermally-conductive ice shell, the bottom heat flux is simply given by the static heat equation in spherical geometry. If, by contrast, convection operates within the ice shell, the bottom heat flux is deduced from scaling laws derived from simulations of mix-heated stagnant-lid convection (Deschamps and Vilella, 2021). In practice, we compute both the conductive and convective heat fluxes, $\phi_{cond}$ and $\phi_{conv}$, and assume that convection operates if $\phi_{conv} \ge \phi_{cond}$.

To relate the ice shell Rayleigh number, $Ra$, which measures the vigor of convection, and the viscous temperature scale $\gamma$, which controls the variation of viscosity with temperature, to the bottom heat flux, we use the following relations

$$ Ra = \frac{\alpha_I \rho_I g \Delta T D^3}{\eta_m \kappa_I}, \quad (4) $$

$$ \gamma = \frac{E^* \Delta T}{R_{gas} T_m^2}, \quad (5) $$
$$ \phi_{bot} = 1.46 \frac{k_I \Delta T}{f^2} \frac{Ra^{0.27}}{\gamma^{1.21}}, \qquad (6) $$

where $\alpha_I$, $\rho_I$, $E^*$, $\kappa_I$, and $k_I$ are the thermal expansion coefficient, density, activation energy, thermal diffusivity, and thermal conductivity of solid ice, respectively, g is the gravitational acceleration of the body, $T_m$ is the interior temperature of the ice shell, $\eta_m$ the ice viscosity at this temperature, D is the ice shell thickness, $\Delta T = T_{bot} - T_{surf}$ the temperature jump across it, f the ratio between its inner and outer radii ($r_{bot}/R$), and $R_{gas}$ is the ideal gas constant. Surface temperature is set to 40 K after McKinnon et al. (1997) and kept constant during the simulations. Numerical values for all parameters are listed in Table 3.

The interior temperature $T_m$ is defined as the averaged temperature within the well-mixed convective interior, i.e., excluding the stagnant lid and thermal boundary layers, and is given by the following relationship (see Deschamps and Vilella, 2021, for details),

$$ T_m = T_{bot} - 1.23 \frac{R_{gas} T_m^2}{E^* f^{1.5}} + (3.5 - 2.3f) \frac{1 + f + f^2}{3} \frac{\rho_I H D^2}{k_I \Delta T} \frac{\Delta T}{Ra^{0.25}}, \quad (7) $$

where H is, as before, the internal heating due to tidal dissipation. Because $Ra$ implicitly depends on $T_m$ (through the viscosity), Eq. (7) does not have an analytical solution. We solved Eq. (7) using a standard Newton-Raphson method. The interior viscosity is deduced from $T_m$ following

$$ \eta_m = \eta_{ref} \exp \left[ \frac{E^*}{R_{gas}T_{ref}} \left( \frac{T_{ref}}{T_m} - 1 \right) \right], \quad (8) $$

where $\eta_{ref}$ is the viscosity of pure water ice close to the melting point, which, in our case, is the temperature of the water liquidus at the bottom of the ice shell $T_{ref} = T_{bottom}$. For $E^* = 60$ kJ/mol, typical of ice Ih (Durham et al., 1997), Eq. (8) generates strong viscosity changes with temperature leading to the formation of a stagnant lid in the top part of the ice shell. The thermal conductivity of ice I, $k_I$, depends on temperature (Anderson and Suga, 1994), but is here assumed constant throughout the ice shell. Deschamps (2021) showed that temperature-dependent conductivity has a strong impact on the thickness and thermal structure of the stagnant lid, but that the thermal evolution of the ice shell is nevertheless well-approximated by a homogeneous conductivity provided that the conductivity is considered at the temperature at the bottom of the ice shell, which, for most icy bodies spans the range 2.0-3.0 W/m/K. Here we assumed a value of 2.6 W/m/K (Deschamps, 2021).

We may also point out that the heat flux given by Eq. (6) does not explicitly depend on H. It is, however, affected by H through $Ra$ and $\gamma$, which both depend on $T_m$ and thus on H. It is also worth noting that, while Eqs. (6) and (7) were built from simulations assuming a homogeneous distribution of H throughout the ice shell, they also describe simulations with viscosity-dependent heating, provided that the maximum tidal dissipation (the highest H) occurs in the hottest regions (Deschamps and Vilella, 2021).

In addition to the amount of internal heating within the ice shell, the presence of impurities within the subsurface ocean can affect the thermal evolution of these bodies (e.g., Deschamps and Sotin, 2001; Grasset and Pargamin, 2005, and references therein) by lowering the melting temperature of the ice shell and thereby changing its crystallization behavior. This is associated with an increase in bulk viscosity, which reduces the vigor of convection and the efficiency of heat transfer, further delaying the crystallization of the ice shell. New Horizons observations of the surfaces of Pluto and Charon reveal CH$_4$, N$_2$, and CO ices (Grundy et al., 2016). In addition, H$_2$O and NH$_3$ ices have also been detected that could have originated in the interior (Dalle Ore et al., 2019). Here, we consider ammonia (NH$_3$) as the main impurity species (e.g., Desch et al., 2009; Waite Jr. et al., 2009; Clark et al., 2014; Grundy et al., 2016; Nimmo and Pappalardo, 2016; Dalle Ore et al., 2018; Beyer et al., 2019) and rely on the liquidus of the H$_2$O-NH$_3$ system of Deschamps and Sotin (2001). Other impurity compounds, in particular magnesium

sulfates and methane, may also be present in icy world objects with qualitatively similar effects as ammonia (Vance et al., 2018; Kamata et al., 2019; Vilella et al., 2020). One thing to bear in mind in the context of the H$_2$O-NH$_3$ system, is that only water crystallizes until the eutectic composition is reached (32.2 vol% NH$_3$). Impurities (NH$_3$) remain in the ocean, whereby their concentrations increase and, as a result, impact the ice shell as it thickens (by increasing viscosity).

Finally, to compute the thermal evolution of Pluto and Charon, we integrate Eq. (1) using an adaptive step-size-controlled 4th-order Runge-Kutta method, and we deduce the thermal properties of the ice shell (including its average temperature and the thickness of the stagnant lid) according to scaling laws developed in Deschamps and Vilella (2021), to which the reader is referred for more details. As for the timing of the formation of the Pluto system, we assume (1) that the accretion process has been completed after the CAIs have formed, and therefore disregard the contribution from short-lived radionuclides, and (2) that Charon formed 200 Myr after the formation of Pluto (e.g., Canup et al., 2021). In connection with the initial thermal state of Pluto and Charon, we assume an initially homogeneous interior temperature, and consider both a cold- and a hot-start case (e.g., Bierson et al., 2020; Renaud et al., 2021). Cold-start (180 K) corresponds to an absent or very thin ocean layer (15 km thickness), whereas a hot-start commences right at the point where the ocean layer (280 km thickness) starts to crystallize, corresponding to ~270 K for pure water and ~250 K in the case contaminants are present. The tidal heating will be discussed in the next section.

## 3.2. Tidal evolution of a highly eccentric non-synchronous rotating binary system

### 3.2.1. Tidal response

An orbiting moon will raise a tide on a planet, resulting in an imposed potential $\Phi$ that causes the planet to deform. The deformation induces a gravitational potential $\psi$ given by

$$ \psi_I(\mathbf{r}) = \left(\frac{R}{r}\right)^{l+1} k_I \Phi_I(\mathbf{R}, \mathbf{r}^*), \quad (9) $$

where R is the radius of the planet, R is a point on the planet's surface, r is a point exterior to R, while r* describes the position of the perturbing body, and both $\Phi$ and $\psi$ are expanded in terms of spherical harmonics of degree l. The $k_I$ are tidal Love numbers of degree l and depend on the interior properties of the deformed body and determine the amplitude of its response (e.g., Efroimsky, 2012).

The expression for the potential of the tidally deformed planet (Eq. (9)) is valid for a purely elastic response of the planet. In this case, the tidal bulge raised on the planet by the moon is aligned with the direction from the planet's center toward the moon, with no lagging. This ensures that the torques, applied by the moon on the planet and the opposite torque with which the planet is acting on the moon, are zero. Consequently, in the elastic case, there is no influence of the planetary tides on the moon's orbital parameters (semi-major axis, eccentricity, and inclination), and, therefore, tidal heat is not generated in the planet.

Realistic objects, however, deviate from the purely elastic case, as a result of which the tidal bulge acquires a complex structure and is no longer aligned with the perturbing body. Following Efroimsky and Makarov (2014), we decompose the bulge over the tidal Fourier modes. This results in harmonics, of which some lag and some are in advance of the sub-satellite point, but each harmonic produces tidal heat. In this case Eq. (9) needs to be modified, in the time domain, to yield

$$ \psi_I(\mathbf{r}, t) = \left(\frac{R}{r}\right)^{l+1} \hat{k}_I \Phi_I(\mathbf{R}, \mathbf{r}^*), \quad (10) $$

where $\hat{k}_I$ is the (linear) Love operator that maps the entire history of the perturbation ($\Phi_I(t')$ over $t' \le t$) on the value of $\psi$ at the present t. In the time domain, this can be written in the form of a convolution, while
in the frequency domain it is given as a product of Fourier components by

$$ \bar{\psi}_l(\mathbf{r}, \omega_{\mathrm{pq}}^{\mathrm{nm}}) = \left(\frac{R}{r}\right)^{l+1} \bar{k}_l(\omega_{\mathrm{pq}}^{\mathrm{nm}}) \bar{\Phi}_l(\mathbf{R}, \mathbf{r}^*, \omega_{\mathrm{pq}}^{\mathrm{nm}}), \quad (11) $$

where $\omega_{\mathrm{pq}}^{\mathrm{nm}}$ are the Fourier tidal modes (whose absolute values are the physical forcing frequencies exerted in the material) and $\{\omega_{\mathrm{pq}}^{\mathrm{nm}}\}$ are integers that are used to number the modes. In the latter expression, overbars are employed to denote complex Fourier components, i.e.,

$$ \bar{k}_l(\omega_{\mathrm{pq}}^{\mathrm{nm}}) = \mathrm{Re}[\bar{k}_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})] + i \mathrm{Im}[\bar{k}_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})] = |\bar{k}_l| e^{-\epsilon_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})}. \quad (12) $$

In Eq. (11), $\bar{\psi}_l(\mathbf{r}, \omega_{\mathrm{pq}}^{\mathrm{nm}})$ is lagging behind $\bar{\Phi}_l(\mathbf{R}, \mathbf{r}^*, \omega_{\mathrm{pq}}^{\mathrm{nm}})$ by the phase angle $\epsilon_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})$, which by convention, is the negative argument of the complex Love number $k_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})$. Replacing $\mathbf{r}$ by the position of the perturbing body $\mathbf{r}^*$, the additional potential “felt” by the latter (Charon) is obtained.

As shown in Efroimsky and Makarov (e.g., 2014), to each tidal mode $\omega_{\mathrm{pq}}^{\mathrm{nm}}$ corresponds an appropriate Fourier contribution that is proportional to the sine of the phase lag at that mode. By convention, the quantity inverse to the absolute value of this sine is termed the tidal quality factor, $Q_l$, and defined as

$$ \frac{1}{Q_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})} = \sin |\epsilon_l(\omega_{\mathrm{pq}}^{\mathrm{nm}})|. \quad (13) $$

In the following, we concentrate on the tides related to degrees 2 and 3. While both $k_l$ and $Q_l$ depend on interior properties (density and rigidity), $Q_l$ is highly sensitive to viscosity and, therefore, to temperature (Eq. (8)). We will briefly describe this in the next section. For the model specified in Table 2, we obtain present-day degree-2 and -3 tidal Love numbers of $k_2 = 0.09$, $h_2 = 0.28$, and $l_2 = 0.058$ and $k_3 = 0.05$, $h_3 = 0.2$, and $l_3 = 0.06$, respectively. Equivalent Love numbers for Charon are $k_2 = 0.001$, $h_2 = 0.002$, and $l_2 = 0.001$ and $k_3 = 0.0006$, $h_3 = 0.0014$, and $l_3 = 0.0004$, respectively. The considerable difference between the Love numbers of the two objects is partly due to the larger size of Pluto and partly related to the presence of a large subsurface ocean within Pluto, which causes Pluto to be more “flexible” than solid Charon.

### 3.2.2. Viscoelastic dissipation

Tidal heating takes place in the planetary bodies as a result of a variety of viscoelastic creep processes (e.g., Jackson and Faul, 2010). Several rheological models (e.g., Maxwell, Andrade, Burgers, and Sundberg-Cooper) have been proposed to model viscoelastic creep based on laboratory measurements (e.g., Jackson and Faul, 2010; Sundberg and Cooper, 2010) and have been employed to model dissipation in planetary bodies (e.g., Roberts and Nimmo, 2008; Harada et al., 2014; Efroimsky, 2015; Williams and Boggs, 2015; McCarthy and Cooper, 2016; Renaud and Henning, 2018; Khan et al., 2018; Bagheri et al., 2019; Tobie et al., 2019). The Sundberg-Cooper viscoelastic model (Sundberg and Cooper, 2010) is used in this study. This model is a composite viscoelastic model that includes features from both of the commonly-used Andrade and extended Burgers models, i.e., the “response broadening” behavior and the experimentally-observed secondary dissipation peak (Jackson and Faul, 2010). The Sundberg-Cooper model has been shown to cover the viscoelastic properties of both ice and silicate materials (Sundberg and Cooper, 2010; Caswell et al., 2015; Caswell and Cooper, 2016; McCarthy and Cooper, 2016). However, laboratory measurements of torsional forced oscillations of silicate materials show that below temperatures of about 900 K, viscoelastic creep does not occur (Jackson and Faul, 2010). In the case of Pluto and Charon, temperatures of the silicate core do not exceed 800 K during the orbital evolution, even in a hot-start formation scenario (Bierson et al., 2020). Consequently, dissipation in the silicate part is negligible with quality factors generally exceeding 500. In contrast, ice is considerably more dissipative than silicate with quality factors as low as 0.1 (McCarthy and Cooper, 2016). The Sundberg-Cooper viscoelastic model is described in more detail in Appendix A.

Finally, tidal dissipation in the liquid ocean layer is assumed negligible. This will be discussed further in Section 4.3. In summary, tidal dissipation in Pluto and Charon is predominantly taking place in their ice shells.

An important parameter of any viscoelastic model is the frequency exponent $\alpha$, which determines the variation of the tidal response with frequency. Higher $\alpha$ corresponds to larger dissipation as the period increases ($Q \sim \omega^{-\alpha}$), implying, for example, that as the two objects recede from one another, the stable synchronous state may be reached in a shorter time interval. The exact value of $\alpha$ is not well-constrained, but most laboratory measurements suggest values in the range 0.25-0.33, which is valid for silicate and icy materials (Jackson and Faul, 2010; McCarthy and Cooper, 2016). Here, we use $\alpha \approx 0.27$, but we will also consider the effect of variations in $\alpha$ on the orbital evolution.

### 3.2.3. Tidal evolution model

After the putative Charon-forming giant impact (Canup, 2005; Stern et al., 2006; Arakawa et al., 2019), the orbits of Pluto and Charon have been separated because of the transfer of angular momentum in analogy with the Earth-Moon system. The tidal evolution of Pluto and Charon consists of expansion of the post-impact orbit of the satellite around the dwarf planet, concomitantly with damping of the initially highly eccentric orbit and despinning of both objects from initially higher spin rates (Table 1).

Here, we extend the tidal evolution model of Bagheri et al. (2021) to include the case of non-synchronous rotation, which builds upon the Darwin-Kaula (Kaula, 1964) and Boué and Efroimsky (2019) tidal models through the use of higher-order eccentricity functions and harmonic modes. As the initial orbit of Charon around Pluto is believed to have been highly eccentric (Canup, 2005), the use of the extended tidal model is essential. Because of the present-day near-spherical shape of both bodies and lack of observed libration, we neglect the contribution from librational tides. Regardless, libration would only act to shorten the “tidal evolution lifetime” relative to the “radiogenic heating lifetime”. Moreover, given the present-day near-zero eccentricity, libration is unimportant.

For a non-synchronous planet hosting a non-synchronous satellite without libration, the tidal rates of the semi-major axis $a$, eccentricity $e$, spin rate $\dot{\theta}$, inclination $i$, and tidally-dissipated energy $E$ can be written in terms of the mean motion ($n$), planet and satellite masses ($M$ and $M'$), the planet and satellite radii ($R$ and $R'$), the quality functions ($K_l = k_l/Q_l$ and $K'_l = k'_l/Q'_l$), and spin rates ($\dot{\theta}$, and $\dot{\theta}'$):

$$ \frac{da}{dt} = \frac{G(M + M')}{n a^2} \frac{M'}{M} \left(\frac{R}{a}\right)^5 \left[ \underbrace{F(K_l, \dot{\theta}, n, e)}_{\text{tides on primary body}} + \underbrace{F(K'_l, \dot{\theta}', n, e)}_{\text{tides on satellite body}} \right], \quad (14) $$

$$ \frac{de}{dt} = \frac{G(M + M')}{n a^3} \frac{M'}{M} \left(\frac{R}{a}\right)^5 \left[ \underbrace{\mathcal{L}(K_l, \dot{\theta}, n, e)}_{\text{tides on primary}} + \underbrace{\mathcal{L}(K'_l, \dot{\theta}', n, e)}_{\text{tides on satellite}} \right], \quad (15) $$

$$ \left(\frac{d\dot{\theta}}{dt}\right)_{\text{primary}} = \frac{GM'^2}{a^3 MR^2} \left(\frac{R}{a}\right)^3 \left[ \underbrace{\mathcal{G}(K_l, \dot{\theta}, n, e)}_{\text{tides on primary}} \right], \quad (16) $$

$$ \left(\frac{di}{dt}\right)_{\text{primary}} = n \sin i \frac{M'}{M} \left(\frac{R}{a}\right)^5 \underbrace{\left[ I(K_l, \dot{\theta}, n, i, e) \right]}_{\text{tides on primary}}, \quad (17) $$

$$ \left(\frac{dE}{dt}\right)_{\text{primary}} = \frac{GM'^2}{a} \left(\frac{R}{a}\right)^5 \underbrace{\left[ \mathcal{H}(K_l, \dot{\theta}, n, e) \right]}_{\text{tides on primary}}, \quad (18) $$

Identical expressions for the time derivatives of $\dot{\theta}$, $i$, and $E$ of the secondary are obtained by replacing $R$, $K_l$, and $\dot{\theta}$ with those of the secondary, i.e., $R'$, $K'_l$, and $\dot{\theta}'$, and interchanging $M$ and $M'$. Detailed expressions for the functions $F$, $L$, $I$, $G$, and $H$ are given in Appendix B.
As apparent from Eqs. (14)–(18), our tidal model includes dissipation within both the primary (Pluto) and the satellite (Charon) contrary to previous studies that only accounted for dissipation in Pluto (Robuchon and Nimmo, 2011; Barr and Collins, 2015; Hammond et al., 2016; Arakawa et al., 2019).

Our tidal model is also favored over the more commonly-used Constant Time Lag (CTL) and Constant Phase Lag (CPL) models (e.g., Heller et al., 2011; Cheng et al., 2014; Samuel et al., 2019; Arakawa et al., 2019; Ferraz-Mello et al., 2020). The CTL model erroneously implies that all the tidal strain modes experience the same temporal delay relative to the modes comprising the tidal stress (Efroimsky and Makarov, 2013; Makarov and Efroimsky, 2013). The CPL model, on the other hand, is not supported by physical principles because it assumes a constant tidal response independent of the rotation frequency, which is not supported by laboratory and geophysical observations (Jackson and Faul, 2010; Jackson, 2005; Khan et al., 2018; Bagheri et al., 2019; Lau and Faul, 2019). The model employed here assigns a separate phase lag and amplitude response with appropriate adjustment for frequency dependence of the tidal mode through the viscoelastic model (Sundberg-Cooper), as discussed above.

Since the Pluto-Charon system is presently locked in a 1:1 spin-orbit resonance, backward integration of the orbit is not possible. Instead, several forward-in-time model computations are conducted starting from different initial conditions and searching for those model runs that lead to the observed orbital properties (Table 1). The heat arising from tidal dissipation is computed using Eq. (18), while the time evolution of the orbital parameters are computed via Eqs. (14)–(17). Finally, the initial inclination of Charon’s orbit is assumed to be zero from the time of formation, which is a typical outcome of impact formation simulations (e.g., Ida et al., 1997; Canup, 2005; Canup and Salmon, 2018; Citron et al., 2015), as a result of which we do not track evolution in inclination.

## 4. Results and discussion

A “nominal case” is chosen among the simulations that matches the present-day observed orbital parameters (Table 1) with interior properties and thermal parameters as given in Tables 2 and 3, respectively. The initial orbital eccentricity, semi-major axis, and spin rates are based on the Pluto-Charon impact simulations of Canup (2005) and reflect the fact that in the post-collisional state, Charon’s orbital eccentricity is high and both bodies commenced with spin rates that are higher than their initial orbital mean motion, indicative of a closer-in satellite. Values employed are $e \sim 0.4$, $a \sim 0.65a_p$, and $\sim 5$ times the initial mean motion, corresponding to a rotation period of $\sim 10$ hr, where $a_p$ is the present-day observed semi-major axis. The choice of initial parameters will be discussed further in Section 4.3.

### 4.1. Tidal evolution

Fig. 2 shows the evolution of the spin rates of the two bodies, semi-major axis of the orbit, orbital eccentricity, and the tidal heat generated in each of the objects. Fig. 2a shows that the tidal evolution of the binary system occurs very rapidly with the stable 1:1 spin-orbit resonance state ($\theta/n = 1$) achieved in $\sim 2 \times 10^5$ years. Moreover, Charon’s spin rate monotonically decreases until it falls into the 3:2 spin-orbit resonance ($\theta/n = 1.5$) after $\sim 10^4$ years. Once the eccentricity of the orbit has become sufficiently damped ($e < 0.2$), the moon rapidly falls into the stable 1:1 spin-orbit resonance. A consequence of the capture into these resonances is a slower orbital evolution of Charon and despinning rate of Pluto. Because of its smaller size, Charon reaches the 1:1 tidally-locked state much faster than Pluto. Pluto’s spin rate, in contrast, increases in the beginning of the evolution but reverts once the satellite reaches the 3:2 resonance and only starts to properly decrease toward synchronization after Charon has entered the 1:1 resonance. The initial spin-up of Pluto can be understood from the point of

view of conservation of angular momentum. As Pluto spins up, Charon’s angular momentum decreases. Since Charon’s orbit is expanding, an eccentricity increase is required. This eccentricity increase, seen in Fig. 2c at around $10^3$ yr, is balanced by a fractional expansion in semi-major axis (Fig. 2b).

Fig. 2b shows the variation in orbital separation. As illustrated in the plot, the orbit does not expand monotonically because of the aforementioned resonances. In fact, and as explained above, the orbit expands for a short period of time just before Charon’s capture into the 3:2 spin-orbit resonance at which point the orbital separation decreases. From the point of view of conservation of angular momentum, the fractional decrease in the observed orbital separation comes about because of the decrease in eccentricity and spin-down of Pluto. Once the satellite falls out of this resonance and has reached the 1:1 spin-orbit resonance, Pluto is propelled toward synchronization. At this point the orbit begins to expand as angular momentum is transferred from Pluto’s super-synchronous spin rate into their mutual orbit, increasing the semi-major axis until the orbital evolution stabilizes at $a = a_p$ when the system has reached its final present-day dual-synchronous state. A qualitatively similar behavior of the orbital evolution of the Pluto-Charon system was observed by Renaud et al. (2021), but because of a different initial $\dot{\theta}/n$ and differences in physical properties of the planetary objects, the exact timing of the 3:2 and 1:1 resonance captures vary slightly in comparison to our results.

Finally, Fig. 2d displays the heat generated in the two bodies due to both tidal and radiogenic heating. We can make a number of observations from this plot: (1) the amount of heat generated by tidal dissipation is higher in Charon than in Pluto, reaching $10^{-7}$ W/kg for Charon and $10^{-8}$ W/kg for Pluto, respectively, and is, relative to radiogenic heating, only of importance in the early stages of orbital evolution ($< 10^3$ yr); (2) radiogenic heating approximately equals the maximum tidal heat generated in Pluto, and is approximately 10 times less in the case of Charon; (3) the tidal heat source rapidly decreases as the eccentricity of the orbit diminishes; (4) radiogenic heating appears to govern the thermal evolution throughout most of the history of both bodies; and (5) as Pluto and Charon reach the dual-synchronous state, there is no longer any tidal force operating and tidal dissipation ceases to be an active heat source. As the subsurface ocean evolves with time, the thermal evidence of despinning is no longer visible. Nonetheless, the tidal dissipation results in stresses within the body (not studied here) that are potentially visible as surface features. If such attributes are able to survive until the present time, they can, in principle, be used to constrain the past history of the objects (e.g., Rhoden et al., 2010; Patthoff and Kattenhorn, 2011; Robuchon and Nimmo, 2011). However, Rhoden et al. (2020) observed a lack of correlation between the tides and surface fractures on Charon. This implies that ocean solidification postdates the synchronization of the orbit, in agreement with the observation made here, that any tidally-induced signature on the surface has been removed.

### 4.2. Evolution of thermal structure

The thermal evolution of Pluto and Charon for the nominal case is shown in Fig. 3. Both nominal cases start with a conservative 15 km ocean thickness (other relevant parameters are compiled in Table 3), and develop oceans (see insets in Fig. 3a,b). In our computations, we focused on the time evolution of the ocean thickness and we found that convection facilitates growth of the ice layer because of the efficient removal of heat, which promotes the crystallization of the ocean. Our model runs actually show that while in most cases oceans develop on both bodies, these are unlikely to persist on Charon until the present. The exact conditions for this to occur depend on several parameters (mainly reference viscosity, radioactive element content, and impurity content) that will be discussed further in Section 4.3.

More specifically, Fig. 3 illustrates the thermal evolution of Pluto (Fig. 3a) and Charon (Fig. 3b) with time in terms of snapshots of
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Tidal evolution of the Pluto–Charon system. (a) spin rates ($\dot{\theta}/n$), (b) semi-major axis ($a/a_p$), (c) eccentricity ($e$), and (d) generated heat (E) per mass (m) from tidal dissipation and decay of long-lived radioactive isotopes. *n* and *a*<sub>p</sub> are the mean anomaly and the present-day distance between Pluto and Charon (see Table 1). 3:2 SOR and 1:1 SOR in panels (a) and (b) refer to the 3:2 and 1:1 spin–orbit-resonances between Pluto and Charon.</figcaption></figure>

the temperature profile (every 200 Myr). The central part made up of the silicate core heats up as a result of radioactive decay reaching temperatures of ~1300 K (Pluto) and ~800 K (Charon) after 4.5 Gyr of evolution. As core temperatures continue to rise, heat escapes and warms up the surrounding ice shell that softens, allowing for increased dissipation. Note that this is only really of importance in the very earliest stages of evolution, since tidal dissipation only dominates in the first $10^4$ yr (Fig. 3d). Oceans grow progressively with time on both bodies as a consequence of the flow of heat from below and the accompanied melting of the ice layer. While the ocean on Pluto reaches a thickness of 100 km, it re-freezes entirely on Charon (Fig. 3c) on account of its size and therefore limited radiogenic budget. As illustrated in Fig. 3a (inset), convection on Pluto starts operating in the ice shell after ~300 Myr, resulting in a thermal structure composed of a hot thermal boundary layer (TBL) at the bottom, an adiabatic region in the middle of the shell, and a cold TBL and stagnant lid at the top, as also observed elsewhere (e.g., McKinnon, 2006; Deschamps et al., 2010; Robuchon and Nimmo, 2011). In comparison, because of Charon's smaller size and thinner ice shell, heat is removed effectively enough through conduction, and convection in the bottom part of the ice layer never sets in, in agreement with previous work (e.g., Bierson et al., 2018; Nimmo et al., 2017). In this particular simulation, Charon's subsurface ocean solidified completely ~300 Myr ago, consistent with previous studies (e.g., Bierson et al., 2018; Rhoden et al., 2020) and the observation of extensional surface tectonic features (see Section 2). Up until about 3.5 Gyr on Charon, only H<sub>2</sub>O crystallizes out of the ocean. As a consequence, the concentration of NH<sub>3</sub> in the ocean grows, reaching the eutectic at ~3.5 Gyr, at which point we observe a change in slope in the crystallization behavior of the ocean, which is due to the fact that both NH<sub>3</sub> and H<sub>2</sub>O start to crystallize and the concentration of NH<sub>3</sub> in the remaining ocean no longer changes. An interesting consequence of this result is that the bottom of the ice shell on Charon may be composed of a mix of H<sub>2</sub>O and NH<sub>3</sub>, assuming that the latter compound was present in the initial ocean. In contrast, on Pluto the eutectic is never reached.

Fig. 3d displays the evolution of surface and core heat flux (evaluated at the surface) on both bodies. As the central parts start to heat up, core heat flux increases, which continuously melts the initially

conductive ice layer. Yet, as the temperature of the core increases, the viscosity of the ice shell decreases, which results in an increase of the convective heat flux (Eq. (8)). A maximum amount of melting is reached after about 1 Gyr with ocean thicknesses of ~110 km and ~50 km for Pluto and Charon, respectively. At this point core and surface heat flux almost balance and this near-steady-state condition dominates the remainder (over ~3.5 Gyr) of the evolution of Pluto and Charon as core heat flux continuously diminishes because of insufficient radioactive heating, and the ocean layer finally refreezes. The spike in heat flux at 200 Myr is due to the tidal heating that is deposited in the ice layer as Pluto despins. This excess energy represents a transient episode, resulting in short-term melting of the ice that re-freezes as the transient fades away.

## 4.3. Influence of model parameters

### 4.3.1. Orbital parameters

In the nominal case (Fig. 2), we relied on the outcome of impact models (e.g., Canup et al., 2021) to guide our initial choice of orbital parameters. In the following, we briefly discuss the impact of these parameters on the tidal energy that is dissipated and the time for the Pluto–Charon system to evolve to the 1:1 synchronous state. Note that, unlike the thermal parameters, the effect of the orbital parameters, i.e., initial *a*, *e*, and $\dot{\theta}$ (hereinafter $a_0$, $e_0$, and $\dot{\theta}_0$), are not simply studied independently of each other. This is due to the fact that by varying one of the parameters, the present-day observed orbital configuration (Table 1) may not be achieved. For example, for an initial eccentricity of 0.2, with $a_0$ and $\dot{\theta}_0$ equal to the nominal case, the final 1:1 spin–orbit resonance state is only reached at $a = 1.2a_p$ (not shown).

In contrast to the nominal case, here we focus on an example, where the initial orbital parameters are closer to the present-day orbital values: $e_0 = 0.06$, $\dot{\theta}_0/n_0 = 3$, and $a_0 = 0.7a_p$ ($n_0$ is the initial mean motion). For this particular case, the maximum tidal heating rate in each of the two bodies is approximately one order of magnitude less than the nominal case, while the total time to reach the stable 1:1 spin–orbit state is not considerably different from that of the nominal case. We also considered the case where the initial spin rates of the two bodies were not equal. While such a case is plausible, we did not
<figure><img src="image_6.png" /><figcaption>Fig. 3. Temporal evolution of the thermal state of Pluto and Charon. (a) Radial temperature profiles of Pluto with each profile representing a time step of 200 Myrs (large panel) and zoom-in of water and ice layers (inset) in time steps of 15 Myrs. Note that in the inset, the convective ice shell starts by increasing (blue–dark green) in thickness but ends up by decreasing toward the end (light green). (b) as in (a) for Charon. Note that no convective ice shell develops on Charon. (c) Temporal evolution of the thickness of the liquid subsurface ocean on Pluto and Charon. (d) Temporal evolution of the surface heat flow on Pluto and Charon. The core heat flow has been scaled to that at the surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

observe any significant difference in the main characteristics of the orbital evolution from the point of view of the time it takes for the Pluto–Charon system, relative to its entire lifetime, to arrive at the final state and in terms of the amount of tidally-generated heat.

As part of our simulations, we also considered cases where $a_0 > a_p$, i.e., initial Planet–satellite separation is greater than their current distance. However, we found that none of them converged to the present-day orbit. This implies that the separation between the host and the satellite initially had to be smaller than the present-day semi-major

axis and that the subsequent orbital evolution has acted to expand the orbit.

We also investigated the effect of the frequency dependence ($\alpha$) on the orbital evolution, which mainly affects the heating rate caused by tidal dissipation. We considered two possibilities: $\alpha = 0.1$ and $\alpha = 0.4$, and found that in both cases the main features of the orbital evolution relative to the nominal case are not significantly different. For $\alpha = 0.4$, the maximum heating rate can increase by up to one order of magnitude. Generally, we find that in all of our simulations the orbital evolution does not exceed $\sim 10^6$ years and, relative to the contribution
from radioactive decay, tidal heating has little effect on the final state of the Pluto-Charon system.

### 4.3.2. Parameters governing the thermal evolution

Our results show that for the nominal case, Pluto likely harbors a present-day ocean overlain by a conductive ice shell, whereas the ocean on Charon, although present up until ~3.5 Gyr, has refrozen completely (Fig. 2c). Both of these results are in agreement with the observations by New Horizons of predominantly extensional and compressional tectonic features on Pluto and Charon, respectively. In the nominal case (Fig. 3), we fixed a number of key parameters that govern the thermal evolution of the two bodies. These include: reference viscosity, core size, initial thermal state, and ocean contaminants. In the following, we briefly discuss the influence of each of these parameters on the evolution of the ocean thickness.

**Reference viscosity.** Convection is very sensitive to viscosity and its temperature variation. A high viscosity opposes the flow, reducing the strength of convection, while temperature-dependent viscosity triggers stagnant-lid convection, which alters the heat transfer. In our computations, the reference viscosity $\eta_{ref}$ controls the bulk viscosity of the ice shell (Eq. (8)). Fig. 4a displays the evolution of ocean thickness on Pluto and Charon for several different values of reference viscosity ($\eta_{ref}$): $10^{12}$, $10^{14}$, and $10^{16}$ Pa s (all other relevant parameters are as in Table 3). Clearly, $\eta_{ref}$ has a strong impact on ocean thickness, in that larger viscosities generally result in a slower cooling as convection weakens, which increases the longevity and the thickness of the ocean layer. For all considered values of $\eta_{ref}$, subsurface oceans on Pluto may survive until the present day, while the initially-formed oceans on Charon have all re-solidified.

For Pluto in the low viscosity case, corresponding to $\eta_{ref} = 10^{12}$ Pa s (yellow line), the ice shell melts until about 0.4 Gyr, leading to a ~70-km thick ocean, which remains roughly constant for the remainder of the evolution. $\eta_{ref} = 10^{14}$ Pa s (orange line) corresponds to the nominal case and has been discussed above. Convection was also observed for $\eta_{ref} = 10^{15}$ Pa s (not shown). However, for $\eta_{ref} = 10^{16}$ Pa s (red line), convection never commences, and melting of the ice sheet is dominant in the first ~1 Gyr. This results in the conductive shell becoming progressively thinner as the core continues to heat up, with the ocean reaching a maximum thickness of ~200 km. With core heat flux subsequently declining as radioactive heating diminishes, the ocean starts to re-freeze and the ice shell thickens. Because the ocean re-freezes at a relatively slow rate, a subsurface ocean on Pluto is able to persist to the present day, which occurs in all of the considered cases. This differs from what was seen in the models of Robuchon and Nimmo (2011), where no ocean was observed to survive until the present when convection sets in, corresponding to reference viscosities <10¹⁵ Pa s. This difference most probably reflects incorporation of impurities (NH₃) in this study, which acts to change the melting point of the ocean (see discussion on impurities below). On Charon, convection is unable to operate, because (1) the ice shell is not thick enough and (2) the effect of temperature-dependent viscosity is too strong, and an ocean, while present for ~3 Gyr, is unable to remain liquid to the present. Due to the absence of convection on Charon the influence of viscosity is considerably smaller than in the case of Pluto.

**Core size and initial ocean/ice layer thickness.** Varying core size implies changing the bulk radioactive element content and, through that, core heat flow. We varied core size on Pluto (820–870 km) and Charon (390–430 km) and adjusted core density (3–3.6 g/cm³) to ensure that the total mass of the body remains unchanged. All other relevant parameters are as in Table 3. The results are shown in Fig. 4b and illustrate, as above, that oceans on Pluto may still exist today, while oceans on Charon have all re-solidified, with the exception of the case of a small (390 km) core that predicts a thin (≤5 km) subsurface ocean on Charon at the present. Moreover, the changes in core size result in subsurface oceans that vary in thickness by ~30 km and ~20 km

on Pluto and Charon, respectively. Again, convection on Charon is generally inhibited by the presence of a smaller hydrosphere.

Observations of Dione by Cassini, which is similar in size to Charon, allows for a lower rocky core density (Zannoni et al., 2020) than considered in our reference model (Table 2). To test how a lower core density affects ocean thickness, we considered a model for Pluto with a rocky core density and radius of 2.5 g/cm³ and 985 km, respectively, which results in a present-day ocean thickness of ~60 km. For Charon, this implies a larger core with a radius of 470 km, and leads to earlier solidification of the ocean (~1.5 Gyr).

We further investigated the effect of initial ocean/ice layer thickness (Fig. 4c), which, as expected, shows the largest differences in the early stages of evolution, where initially thick ice layers (thin ocean layers) follow the path described earlier (Fig. 4a) for the low reference viscosity cases ($\eta_{ref}=10^{12}$ Pa s), whereas thin ice layers (thick ocean layers) immediately start off by freezing until a quasi-steady-state condition is reached, where the heat entering from below is conductively removed through the ice layer. Independently of initial ocean/ice layer thickness, Fig. 4c shows that subsurface oceans on Pluto and Charon reach thicknesses of ~100 km and 0 km, respectively. These results show that thermal evolution completely removes the signature of the initial ice shell. In summary, the initial ocean/ice layer thickness is less important in the context of the long-term evolution of Pluto and Charon.

**Impurity content.** The presence of impurity species reduces the temperature at the bottom of the ice shell, which increases the bulk viscosity and thus reduces the vigor of convection. To evaluate the importance of this effect on the evolution of Pluto and Charon, we vary the initial fraction of ammonia in the range 1–5 wt% (all other relevant parameters are as in Table 3). The results are shown in Fig. 4d for both bodies. Two conclusions may be drawn: (1) regardless of impurity content, a subsurface ocean is always present on Pluto at the present day, and (2) oceans also form on Charon, but completely solidify after ~3.5–4 Gyr. Between the various cases, ocean thicknesses on Pluto and Charon range from ~10–50 km and ~40–120 km, respectively. Convection only happens for $NH_3=1$ and 3 wt%, respectively, whereas for $NH_3=5$ wt%, the ice layer remains conductive. Overall, convection within Charon appears to be difficult to set up and to maintain. For instance, even for $NH_3=1$ wt%, convection is only operative between ~0.25 and 2.4 Gyr. With reference viscosity, impurity content has the largest influence on final ocean thickness, in particular for Pluto, because of its impact on the vigor of convection.

## 5. Summary and outlook

In this work we have implemented a comprehensive semi-analytical tidal model based on the latest advances in tidal theory that considers dual dissipation in a binary system and invokes an appropriate viscoelastic rheology (Sundberg-Cooper) for proper computation of the tidal response functions. Because the spin rate of a planet can impart orbital changes that lead, via dissipation, to changes in the thermal state of the other body, a thermal-orbital feedback exists that calls for a joint approach. Consequently, here we studied the combined tidal and thermal evolution of the Pluto-Charon system with particular emphasis on the possibility for present-day subsurface oceans on both bodies.

We found that a subsurface ocean always develops and remains liquid to the present day on Pluto. Subsurface oceans on Charon also developed in all of the studied cases in its past (up until ~3.5 Gyr), but because of its smaller size and therefore lower radiogenic budget, the ocean re-freezes between 0–500 Myr ago. These observations were found to be robust across an entire range of orbital and thermal model parameter values that were considered in this study. Based on the parameter analysis, the present-day ocean on Pluto ranges in thickness from 40 km to 150 km. These observations are generally consistent with previous modeling results (Robuchon and Nimmo, 2011; Hammond
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Summary of thermal parameter study illustrating their influence on the evolution (ocean thickness) of the Pluto–Charon system. (a) Ice reference viscosity; (b) Core radius; (c) Initial ocean thickness; and (d) Impurity content. In plots (a–d) Pluto is described by yellow, orange, and red lines, whereas Charon is delineated by light gray, gray, and black lines. The change in slope that is observed on all Charon-related curves toward the end of the evolution (~3.5 Gyr) is indicative of a change in the crystallization behavior of the ocean. See main text for more details. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

et al., 2016; Nimmo et al., 2016; Desch and Neveu, 2017; Bierson et al., 2018), but differ in detail mostly because we incorporate impurities. Observations of extensional and compressional surface tectonics on Charon and Pluto, respectively, made by New Horizons support this. While we applied our tidal-thermal model to the Pluto–Charon system, it stands to reason that it can be applied to any satellite orbiting a central body (e.g., Neptune–Triton, Earth–Moon, or an exoplanet system).

When viewed over a time span of ~4.5 Gyr, the tidal energy that is dissipated in both bodies as the orbits evolve toward synchronization, including higher spin–orbit resonances, is less significant relative to the heat emanating from the decay of long-lived radioactive isotopes. Generally, the tidally-deposited energy from despinning is only of relevance in the early stages of the evolution of the Pluto–Charon system. In spite of the fact that tidal heating plays a less important role, our calculations nevertheless show that during the short period that tides were active, the amount of energy released by tides was much higher than that emanating from radioactive decay. This is of particular importance in the case of the tidally-active Jovian and Saturnian satellites (Asphaug and Reufer, 2013; Shoji et al., 2014; Neveu and Rhoden, 2019; Bierson and Steinbrügge, 2021; Lainey et al., 2020; Tyler et al., 2015) and exoplanetary systems (Miller et al., 2009; Shoji and Kurita, 2014; Driscoll and Barnes, 2015; Dobos et al., 2019).

In the context of tidally active and potentially habitable worlds, heating due to tidal dissipation in the oceans, i.e., from friction generated at the ice–water and/or water–solid mantle interfaces, may be important (Tyler, 2008; Beuthe, 2019; Rovira-Navarro et al., 2019). On Earth, for example, the dissipation due to tidally-induced motion in the oceans is greater than that arising from the tidal deformation of the solid-body interior (Egbert and Ray, 2003). In the case of subsurface oceans, resonantly forced tidal states in the ocean with highly elevated power levels may even dominate the viscoelastic tidal dissipation in the solid parts of the body (Tyler, 2020). Tidal dissipation in the liquid ocean is a complex process and depends on parameters such as depth, composition, and the topography of the ice and ocean bottom (Tyler, 2008) that are currently unknown. Here we have excluded the possibility for this additional heat source, but we may observe that had we incorporated the effect of ocean tidal heating, this would act to shorten

the duration for complete synchronization. However, as tidal heating is only really significant very early on, the potential effect of considering ocean tides would not lead to any change in the long-term evolution. From the point of view of tidal heating and evolution timescales, our results represent a lower and an upper bound, respectively.

Because direct information on the interiors of icy ocean worlds is limited, their astrobiological potential is difficult to assess in detail. For this, the key parameters relating to habitability, which include ocean depth, temperature, and chemistry, among others, need to be well-established (e.g. Heller et al., 2021). Given the high resolving power of geophysical methods, particularly of seismology, in situ geophysical exploration of icy planets and satellites is clearly the means to move beyond the current impasse as pointed out in numerous studies (e.g., Kovach and Chyba, 2001; Lee et al., 2003; Panning et al., 2006, 2018; Vance et al., 2018; Stähler et al., 2018; Maguire et al., 2021). The successful Mars Insight mission has paved the way for single-station emplacement and sounding of icy worlds, in that InSight has shown that with a single station, marsquakes can be recorded, located, and used to sound its interior (Knapmeyer-Endrun et al., 2021; Khan et al., 2021; Stähler et al., 2021). In this context, the proposed Europa (Pappalardo et al., 2013) and the selected Titan relocatable lander “Dragonfly” (Barnes et al., 2021) missions both include a single-station seismometer in their proposed scientific payload that should enable better characterization of ocean thickness.

## Acknowledgments

We are grateful to two anonymous reviewers for constructive comments that led to an improved manuscript. This work was supported by a grant from the Swiss National Science Foundation (project 172508). We would like to acknowledge our colleagues Francis Nimmo, Joe Renaud, Michael Efroimsky, and Robert Tyler for informed discussions on aspects of tidal evolution and dissipation.

## Appendix A. Sundberg-Cooper viscoelastic model
<table><caption>**Table 4**
Parameters associated with the Sundberg-Cooper viscoelastic model.</caption><thead><tr><th>Parameter (unit)</th><th>Symbol</th><th>Value/expression</th></tr></thead><tbody><tr><td>Constant</td><td>S</td><td>α! sin(απ/2)</td></tr><tr><td>Constant</td><td>C</td><td>α! cos(απ/2)</td></tr><tr><td>Burgers coefficient</td><td>λ</td><td>(J<sub>U</sub>η<sub>1</sub>)<sup>2</sup> + (J<sub>U</sub>/J<sub>R</sub>)<sup>2</sup></td></tr><tr><td>Relaxed compliance (Pa s)</td><td>J<sub>R</sub></td><td>0.2J<sub>U</sub></td></tr><tr><td>Voigt-Kelvin viscosity (Pa s)</td><td>η<sub>1</sub></td><td>0.02η</td></tr></tbody></table>

The compliance or creep function associated with the Sundberg-Cooper model is (Sundberg and Cooper, 2010):

$$J(\omega) = J_U - \frac{i}{\eta\omega} + \frac{jJ_R}{i - J_R\eta_1\omega} + J_U(iJ_U\eta\omega)^{-\alpha}\alpha!, \quad (19)$$

where $J_U$ and $J_R$ are unrelaxed and relaxed shear moduli, respectively, $\eta_1$ is the Kelvin-Voigt viscosity, $\eta$ is viscosity ($\eta$ for the various layers is listed in Table 3), $\omega$ is frequency, and $\alpha$ is the frequency exponent as defined in the main text (see Section 3.2.2). The frequency-dependent complex shear modulus is given by

$$\Re[G(\omega)] = \frac{1}{J_U H(\omega)} \left[ \left( [J_U \eta \omega]^\alpha C + 1 \right) \lambda + \frac{J_U}{J_R} \right], \quad (20)$$

$$\Im[G(\omega)] = \frac{\eta\omega}{H(\omega)} \left[ (J_U \eta \omega)^{-1-\alpha} S \lambda + \frac{\eta_1}{\eta} + \left( \frac{\eta_1}{\eta} \right)^2 + (\eta \omega J_R)^{-2} \right], \quad (21)$$

where $\lambda$ is the Burgers coefficient, $S$ and $C$ are constants, and the function $H(\omega)$ is defined as

$$H(\omega) = \frac{1}{(J_R \eta \omega)^2} + \lambda + 2 \frac{J_U}{J_R} + \left(\frac{\eta_1}{\eta} + 1\right)^2 + \lambda (J_U \eta \omega)^{-2\alpha} (a!)^2 + 2 (J_U \eta \omega) \left[ S \left( J_U \eta_1 \omega \left(1 + \frac{\eta_1}{\eta}\right) + \frac{J_U}{J_R (J_R \eta \omega)} \right) + C \left(\lambda + \frac{J_U}{J_R}\right) \right]. \quad (22)$$

The parameters and constants used in the above expressions are given in Table 4. The intrinsic shear quality factor $Q^{-1}$, which is not to be confused with the $Q^{-1}$ associated with global tidal dissipation (Eq. (13)), is a measure of dissipation and can be computed from

$$Q^{-1} = \frac{\Im[G(\omega)]}{\Re[G(\omega)]}. \quad (23)$$

## Appendix B. Orbital evolution theory

In what follows, we rely on Bagheri et al. (2021) for $da/dt$ and $de/dt$, and augment these with full expressions for $d\dot{\theta}/dt$, $dE/dt$, and the case of non-synchronous rotation. The time evolution of each of the aforementioned orbital parameters of the two-body system can be cast as

$$\left(\frac{dx}{dt}\right) = \left(\frac{dx}{dt}\right)_{\text{primary}} + \left(\frac{dx}{dt}\right)_{\text{secondary}}, \quad (24)$$

where $x$ is an orbital parameter. The two terms refer to the tides in the planet (primary) and the tides in the satellite (secondary), respectively. Because of the spherical shape of both Pluto and Charon, we can omit the contribution from libration.

The expression for the semi-major axis rate is

$$\frac{da}{dt} = -2an \sum_{l=2}^{\infty} \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \sum_{p=0}^{l} \sum_{q=-\infty}^{\infty} G_{lpq}^{2}(e)(l-2p+q) \left[ \left(\frac{R}{a}\right)^{2l+1} \frac{M'}{M} F_{lm}^{2}(i) K_{l}(\omega) + \left(\frac{R}{a}\right)^{2l+1} \frac{M}{M'} F_{lm}^{2}(i') K_{l'}(\omega') \right], \quad (25)$$

where $M$ and $M'$ are the planet and satellite masses, $R$ and $R'$ are their radii, $K_l = k_l/Q_l$ and $K_{l'} = k_{l'}/Q_{l'}$ are the planet and satellite quality functions, $G(e)$ and $F(i)$ are the eccentricity and inclination functions, respectively, $\omega$ is the tidal mode, $n$ is the mean motion, and all the other variables are as defined in Eqs. (14)–(16). Since the inclinations

of the orbits remain small, only $F_{201}$ and $F_{220}$ are relevant and equal to $\frac{1}{2}$ and 3, respectively. Similarly to $da/dt$, the general expression for the eccentricity rate is

$$\frac{de}{dt} = -\frac{1 - e^2}{e} \frac{n}{M M'} \sum_{l=2}^{\infty} \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \sum_{p=0}^{l} \sum_{q=-\infty}^{\infty} G_{lpq}^2(e)(l-2p+q) \left[ \left(\frac{R}{a}\right)^{2l+1} \frac{M'}{M} F_{lm}^2(i) K_l(\omega) + \left(\frac{R'}{a}\right)^{2l+1} \frac{M}{M'} F_{lm}^2(i') K_{l'}(\omega') \right] \\ - \frac{\sqrt{1 - e^2}}{e} \frac{n}{M M'} \sum_{l=2}^{\infty} \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \sum_{p=0}^{l} \sum_{q=-\infty}^{\infty} G_{lpq}^2(e)(l-2p) \left[ \left(\frac{R}{a}\right)^{2l+1} \frac{M'}{M} F_{lm}^2(i) K_l(\omega) + \left(\frac{R'}{a}\right)^{2l+1} \frac{M}{M'} F_{lm}^2(i') K_{l'}(\omega') \right]. \quad (26)$$

Equations for the tidal heating and the tidal torque are given by

$$\left(\frac{dE}{dt}\right)_{\text{primary}} = \frac{GM'^2}{a} \sum_{l=2}^{\infty} \left(\frac{R}{a}\right)^{2l+1} \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \\ \times \sum_{p=0}^{l} \sum_{q=-\infty}^{\infty} G_{lpq}^2(e) F_{lm}^2(i) K_l(\omega')\omega', \quad (27)$$

$$\left(\frac{d\dot{\theta}}{dt}\right)_{\text{primary}} = \frac{GM'^2}{Ca} \sum_{l=2}^{\infty} \left(\frac{R}{a}\right)^{2l+1} \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \\ \times \sum_{p=0}^{l} \sum_{q=-\infty}^{\infty} G_{lpq}^2(e) F_{lm}^2(i) m K_l(\omega'), \quad (28)$$

where $C$ is the polar moment of inertia of the planet. Similar expressions are obtained for the satellite by replacing $M$ with $M'$ and $C$ by $C'$, where $C'$ is the polar moment of inertia of the secondary (see below).

Due to the slow convergence of the series and the relatively high eccentricities found in this study, we have to include higher-order eccentricity terms to ensure precision of results and stability of integration. In Eq. (24), the contribution of tides raised by the satellite in the planet is

$$\left(\frac{da}{dt}\right)_{\text{primary}} = n \left(\frac{R^5}{a^4}\right) \left(\frac{M'}{M}\right) \times \mathcal{F}(K_l, \dot{\theta}, n, e). \quad (29)$$

The contribution due to the tides raised by the planet in the satellite looks similar

$$\left(\frac{da}{dt}\right)_{\text{secondary}} = n \left(\frac{R'^5}{a^4}\right) \left(\frac{M}{M'}\right) \times \mathcal{F}(K'_l, \dot{\theta}', n, e). \quad (30)$$

Here $\mathcal{F}$ is a function of eccentricity, mean motion, spin rate, and tidal response of either body

$$\mathcal{F}(K_l, \dot{\theta}, n, e) = \sum_{i=0}^{9} e^{2i} \left( \sum_{j=-7}^{-1} K_l \left(jn - 2\dot{\theta}\right) \phi_j^{2i} + \sum_{j=1}^{11} K_l \left(jn - 2\dot{\theta}\right) \phi_j^{2i} + \sum_{j=1}^{9} K_l \left(jn\right) \phi_j^{2i} \right). \quad (31)$$

The coefficients $\phi_j^{2i}$ and $\dot{\phi}_j^{2i}$ of the series are tabulated in Supplementary Tables 1 and 2. Note that in these tables, the terms of the series that are not specifically mentioned are equal to zero. The above equations have been derived for the general case, i.e., with neither of the bodies assumed synchronous. In the specific situation of a synchronized moon, we have $\dot{\theta}' = n$, wherefore the semi-diurnal term in Eq. (30) vanishes: $K_l(2n - 2\dot{\theta}') = 0$. In the contribution from the planet, the semi-diurnal term vanishes when the satellite is at the synchronous orbit, i.e., when $n = \dot{\theta}'$.

Similarly, to compute the eccentricity evolution, we write down all the inputs entering Eq. (26). The input generated by the tides in the planet is

$$\left(\frac{de}{dt}\right)_{\text{primary}} = -n \frac{M'}{M} \left(\frac{R}{a}\right)^5 \times \mathcal{L}(K_l, \dot{\theta}, n, e), \quad (32)$$
while the input from the tides in the satellite is given by

$$ \left(\frac{de}{dt}\right)_{\text{secondary}} = -n \frac{M}{M'} \left(\frac{R'}{a}\right)^5 \times \mathcal{L}(K'_l, \dot{\theta}', n, e), \quad (33) $$

where $\mathcal{L}$ is a function of eccentricity, mean motion, spin rate, and tidal response of either body

$$ \mathcal{L}(K_l, \dot{\theta}, n, e) = \sum_{l=1}^{9} e^{2l-1} \left( \sum_{j=-7}^{-1} K_l (jn - 2\dot{\theta}) \lambda_j^{2l-1} + \sum_{j=1}^{11} K_l (jn - 2\dot{\theta}) \lambda_j^{2l-1} + \sum_{j=1}^{9} K_l (jn) \dot{\lambda}_j^{2l-1} \right). \quad (34) $$

The coefficients $\lambda_j^{2l-1}$ and $\dot{\lambda}_j^{2l-1}$ are tabulated in Supplementary Tables 3 and 4. Note that, similarly to $da/dt$, the expression is general, in that neither of the bodies is a priori assumed synchronous. For the synchronized case, the term associated with the semi-diurnal tide in Eq. (32) vanishes.

To compute the time evolution of the tidally-dissipated energy in the primary body, we expand Eq. (27) as follows

$$ \left(\frac{dE}{dt}\right)_{\text{primary}} = \frac{GM'^2}{a} \left(\frac{R}{a}\right)^5 \times \mathcal{H}(K_l, \dot{\theta}, n, e), \quad (35) $$

whereas the time evolution of the tidally-dissipated energy in the secondary is given by a similar expression

$$ \left(\frac{dE}{dt}\right)_{\text{secondary}} = \frac{GM^2}{a} \left(\frac{R'}{a}\right)^5 \times \mathcal{H}(K'_l, \dot{\theta}', n, e). \quad (36) $$

The expression for $\mathcal{H}$ is a function of eccentricity, mean motion, spin rate, and tidal response of either body

$$ \mathcal{H}(K_l, \dot{\theta}, n, e) = \sum_{l=0}^{9} e^{2l} \left[ \sum_{j=-7}^{-1} \xi_j^{2l} |jn - 2\dot{\theta}| K_l (|jn - 2\dot{\theta}|) + \sum_{j=1}^{11} \xi_j^{2l} |jn - 2\dot{\theta}| K_l (|jn - 2\dot{\theta}|) + \sum_{j=1}^{9} jn \xi_j^{2l} K_l (jn) \right]. \quad (37) $$

The coefficients $\xi_j^{2l}$ and $\dot{\xi}_j^{2l}$ are tabulated in Supplementary Tables 5 and 6.

Finally, the expression for the time evolution of the spin rate of the primary, Eq. (28) can be expanded as

$$ \left(\frac{d\dot{\theta}}{dt}\right)_{\text{primary}} = \frac{GM'}{Ca} \left(\frac{R}{a}\right)^5 \times \mathcal{G}(K_l, \dot{\theta}, n, e), \quad (38) $$

while that for the secondary is obtained from the following expression

$$ \left(\frac{d\dot{\theta}}{dt}\right)_{\text{secondary}} = \frac{GM}{C'a} \left(\frac{R'}{a}\right)^5 \times \mathcal{G}(K'_l, \dot{\theta}', n, e). \quad (39) $$

The expression $\mathcal{G}$ is a function of eccentricity, mean motion, spin rate, and tidal response of either body

$$ \mathcal{G}(K_l, \dot{\theta}, n, e) = \sum_{l=0}^{9} e^{2l} \left[ \sum_{j=-7}^{-1} \chi_j^{2l} K_l (jn - 2\dot{\theta}) + \sum_{j=1}^{11} \chi_j^{2l} K_l (jn - 2\dot{\theta}) \right]. \quad (40) $$

The coefficients $\chi_j^{2l}$ are tabulated in Supplementary Table 7.

## Appendix C. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.icarus.2021.114871.

## References

Anderson, O., Suga, H., 1994. Thermal conductivity of the ih and XI phases of ice. Phys. Rev. B 50, 6583-6588.

Arakawa, S., Hyodo, R., Genda, H., 2019. Early formation of moons around large trans-Neptunian objects via giant impacts. Nat. Astron. 3 (9), 802-807.

Arakawa, S., Hyodo, R., Shoji, D., Genda, H., 2021. Tidal evolution of the eccentric moon around dwarf planet (225088) gonggong. arXiv preprint arXiv:2108.08553.

Asphaug, E., Reufer, A., 2013. Late origin of the Saturn system. Icarus 223 (1), 544-565.

Bagheri, A., Khan, A., Al-Attar, D., Crawford, O., Giardini, D., 2019. Tidal response of Mars constrained from laboratory-based viscoelastic dissipation models and geophysical data. J. Geophys. Res. Planets.

Bagheri, A., Khan, A., Efroimsky, M., Kruglyakov, M., Giardini, D., 2021. Dynamical evidence for phobos and deimos as remnants of a disrupted common progenitor. Nat. Astron. 1-5.

Barnes, J.W., Turtle, E.P., Trainer, M.G., Lorenz, R.D., MacKenzie, S.M., Brinckerman, W.B., Cable, M.L., Ernst, C.M., Freissinet, C., Hand, K.P., Hayes, A.G., Hörst, S.M., Johnson, J.R., Karkoschka, E., Lawrence, D.J., Gall, A.L., Lora, J.M., McKay, C.P., Miller, R.S., Murchie, S.L., Neish, C.D., Newman, C.E., Núñez, J., Panning, M.P., Parsons, A.M., Peplowski, P.N., Quick, L.C., Radebaugh, J., Rafkin, S.C.R., Shiraishi, H., Soderblom, J.M., Sotzen, K.S., Stickle, A.M., Stofan, E.R., Szopa, C., Tokano, T., Wagner, T., Wilson, C., Yingst, R.A., Zacny, K., Stähler, S.C., 2021. Science goals and objectives for the dragonfly titan rotorcraft relocatable lander. Planet. Sci. J. 2 (4), 130.

Barr, A.C., Collins, G.C., 2015. Tectonic activity on Pluto after the Charon-forming impact. Icarus 246, 146-155.

Beuthe, M., 2019. Enceladus's crust as a non-uniform thin shell: II tidal dissipation. Icarus 332, 66-91.

Beyer, R.A., Nimmo, F., McKinnon, W.B., Moore, J.M., Binzel, R.P., Conrad, J.W., Cheng, A., Ennico, K., Lauer, T.R., Olkin, C., Robbins, S., Schenk, P., Singer, K., Spencer, J.R., Stern, S.A., Weaver, H., Young, L., Zangari, A.M., 2017. Charon tectonics. Icarus 287, 161-174.

Beyer, R.A., Spencer, J.R., McKinnon, W.B., Nimmo, F., Beddingfield, C., Grundy, W., Ennico, K., Keane, J.T., Moore, J.M., Olkin, C.B., et al., 2019. The nature and origin of Charon's smooth plains. Icarus 323, 16-32.

Bierson, C., Nimmo, F., McKinnon, W., 2018. Implications of the observed Pluto-Charon density contrast. Icarus 309, 207-219.

Bierson, C.J., Nimmo, F., Stern, S.A., 2020. Evidence for a hot start and early ocean formation on Pluto. Nat. Geosci. 13 (7), 468-472.

Bierson, C.J., Steinbrügge, G., 2021. Tidal heating did not dry out io and Europa. Planet. Sci. J. 2 (3), 89.

Boué, G., Efroimsky, M., 2019. Tidal evolution of the Keplerian elements. Celestial Mech. Dynam. Astronom. 131 (7), 30.

Brozović, M., Showalter, M.R., Jacobson, R.A., Buie, M.W., 2015. The orbits and masses of satellites of Pluto. Icarus 246, 317-329.

Buie, M.W., Tholen, D.J., Grundy, W.M., 2012. The orbit of Charon is circular. Astron. J. 144 (1), 15.

Canup, R.M., 2005. A giant impact origin of Pluto-Charon. Science 307 (5709), 546-550.

Canup, R.M., 2010. On a giant impact origin of Charon, Nix, and Hydra. Astron. J. 141 (2), 35.

Canup, R.M., Kratter, K.M., Neveu, M., 2021. On the origin of the Pluto system. In: The Pluto System After New Horizons. University of Arizona Press, p. 475.

Canup, R., Salmon, J., 2018. Origin of phobos and deimos by the impact of a vesta-to-ceres sized body with mars. Sci. Adv. 4 (4), eaar6887.

Castillo-Rogez, J.C., Lunine, J.I., 2010. Evolution of Titan's rocky core constrained by Cassini observations. Geophys. Res. Lett. 37 (20).

Caswell, T., Cooper, R., 2016. Convection-induced microstructure and tidal dissipation in polycrystalline ice; an experimental approach. In: Lunar And Planetary Science Conference. (1903), p. 1129.

Caswell, T.E., Cooper, R.F., Goldsby, D.L., 2015. The constant-hardness creep compliance of polycrystalline ice. Geophys. Res. Lett. 42 (15), 6261-6268.

Cheng, W., Lee, M.H., Peale, S., 2014. Complete tidal evolution of Pluto-Charon. Icarus 233, 242-258.

Choblet, G., Tobie, G., Sotin, C., Běhounková, M., Čadek, O., Postberg, F., Souček, O., 2017. Powering prolonged hydrothermal activity inside Enceladus. Nat. Astron. 1 (12), 841.

Citron, R.I., Genda, H., Ida, S., 2015. Formation of phobos and deimos via a giant impact. Icarus 252, 334-338.

Clark, R.N., Swayze, G.A., Carlson, R., Grundy, W., Noll, K., 2014. Spectroscopy from space. Rev. Mineral. Geochem. 78 (1), 399-446.

Dalle Ore, C.M., Cruikshank, D.P., Protopapa, S., Scipioni, F., McKinnon, W.B., Cook, J.C., Grundy, W.M., Schmitt, B., Stern, S.A., Moore, J.M., Verbiser, A., Parker, A.H., Singer, K.N., Umurhan, O.M., Weaver, H.A., Olkin, C.B., Young, L.A., Ennico, K., 2019. Detection of ammonia on Pluto's surface in a region of geologically recent tectonism. Sci. Adv. 5 (5).

Dalle Ore, C.M., Protopapa, S., Cook, J., Grundy, W., Cruikshank, D., Verbiser, A., Ennico, K., Olkin, C., Stern, S., Weaver, H., et al., 2018. Ices on Charon: Distribution of H2O and NH3 from new horizons LEISA observations. Icarus 300, 21-32.

Davaille, A., Jaupart, C., 1993. Transient high-Rayleigh-number thermal convection with large viscosity variations. J. Fluid Mech. 253, 141-166.

Desch, S.J., 2015. Density of Charon formed from a disk generated by the impact of partially differentiated bodies. Icarus 246, 37-47.

Desch, S.J., Cook, J.C., Doggett, T., Porter, S.B., 2009. Thermal evolution of Kuiper belt objects, with implications for cryovolcanism. Icarus 202 (2), 694-714.
Desch, S., Neveu, M., 2017. Differentiation and cryovolcanism on charon: A view before and after new horizons. Icarus 287, 175–186.

Deschamps, F., 2021. Stagnant lid convection with temperature-dependent thermal conductivity and the thermal evolution of icy worlds. Geophys. J. Int. 224 (3), 1870–1889.

Deschamps, F., Sotin, C., 2001. Thermal convection in the outer shell of large icy satellites. J. Geophy. Res. Planet. 106 (E3), 5107–5121.

Deschamps, F., Tackley, P.J., Nakagawa, T., 2010. Temperature and heat flux scalings for isoviscous thermal convection in spherical geometry. Geophys. J. Int. 182, 137–154.

Deschamps, F., Vilella, K., 2010. Scaling laws for mixed-heated stagnant lid convection and application to europa. JGR Planets in press.

Dobos, V., Barr, A.C., Kiss, L.L., 2019. Tidal heating and the habitability of the TRAPPIST-1 exoplanets. Astron. Astrophys. 624, A2.

Dobrovolskis, A., Peale, S., Harris, A., 1997. Dynamics of the Pluto-Charon binary. Pluto Charon 159.

Driscoll, P.E., Barnes, R., 2015. Tidal heating of earth-like exoplanets around M stars: thermal, magnetic, and orbital evolutions. Astrobiology 15 (9), 739–760.

Durham, W.B., Kirby, S.H., Stern, L.A., 1997. Creep of water ice at planetary conditions: A compilation. J. Geophy. Res. Planet. 102 (E7), 16293–16302.

Durham, W.B., Prieto-Ballesteros, O., Goldsby, D.L., Kargel, J.S., 2010. Rheological and thermal properties of icy minerals. Space Sci. Rev. 153, 273–298.

Efraimsky, M., 2012. Tidal dissipation compared to seismic dissipation: In small bodies, earths, and super- earths. Astrophys. J. 746, 150.

Efraimsky, M., 2015. Tidal evolution of asteroidal binaries. Ruled by viscosity. Ignorant of rigidity. Astron. J. 150 (4), 98.

Efraimsky, M., Makarov, V.V., 2013. Tidal friction and tidal lagging. applicability limitations of a popular formula for the tidal torque. Astron. J. 764, 26.

Efraimsky, M., Makarov, V.V., 2014. Tidal dissipation in a homogeneous spherical body. I. Methods. Astron. J. 795 (1), 6.

Egbert, G.D., Ray, R.D., 2003. Semi-diurnal and diurnal tidal dissipation from TOPEX/Poseidon altimetry. Geophys. Res. Lett. 30 (17).

Feistel, R., Wagner, W., 2006. A new equation of state for H 2 O ice Ih. J. Phys. Chem Ref. Data 35 (2), 1021–1047.

Ferraz-Mello, S., Beaugé, C., Folonier, H.A., Gomes, G.O., 2020. Tidal friction in satellites and planets. The new version of the creep tide theory. Eur. Phys. J. Spec. Top. 229, 1441–1462.

Grasset, O., Pargamin, J., 2005. The ammonia-water system at high pressures: Implications for the methane of Titan. Planet. Space Sci. 53 (4), 371–384.

Grasset, O., Sotin, C., 1996. The cooling rate of a liquid shell in Titan's interior. Icarus 123 (1), 101–112.

Grundy, W., Binzel, R., Buratti, B., Cook, J., Cruikshank, D., Dalle Ore, C., Earle, A., Ennico, K., Howett, C., Lunsford, A., et al., 2016. Surface compositions across Pluto and Charon. Science 351 (6279).

Hammond, N.P., Barr, A.C., Parmentier, E.M., 2016. Recent tectonic activity on Pluto driven by phase changes in the ice shell. Geophys. Res. Lett. 43 (13), 6775–6782.

Harada, Y., Goossens, S., Matsumoto, K., Yan, J., Ping, J., Noda, H., Haruyama, J., 2014. Strong tidal heating in an ultralow-viscosity zone at the core-mantle boundary of the Moon. Nat. Geosci. 7, 569–572.

Heller, R., Duda, J.-P., Winkler, M., Reitner, J., Gizon, L., 2021. Habitationality of the early earth: liquid water under a faint young sun facilitated by strong tidal heating due to a closer Moon. PalZ 1–13.

Heller, R., Leconte, J., Barnes, R., 2011. Tidal obliquity evolution of potentially habitable planets. Astron. Astrophys. 528, A27.

Hussmann, H., Choblet, G., Lainey, V., Matson, D.L., Sotin, C., Tobie, G., Van Hoolst, T., 2010. Implications of rotation, orbital states, energy sources, and heat transport for internal processes in icy satellites. Space Sci. Rev. 153 (1), 317–348.

Hussmann, H., Sohl, F., Spohn, T., 2006. Subsurface oceans and deep interiors of medium-sized outer planet satellites and large trans-neptunian objects. Icarus 185 (1), 258–273.

Hussmann, H., Sotin, C., Lunine, J.I., 2007. Interiors and evolution of icy satellites, in planets and moons. Treatise Geophys. 10, 509–539.

Hussmann, H., Spohn, T., 2004. Thermal-orbital evolution of Io and Europa. Icarus 171 (2), 391–410.

Ida, S., Canup, R.M., Stewart, G.R., 1997. Lunar accretion from an impact-generated disk. Nature 389 (6649), 353–357.

Jackson, I., 2005. Laboratory measurement of seismic wave dispersion and attenuation at high pressure and temperature. In: Advances In High-Pressure Technology For Geophysical Applications. Elsevier, pp. 95–119.

Jackson, I., Faul, U.H., 2010. Grainsize-sensitive viscoelastic relaxation in olivine: Towards a robust laboratory-based model for seismological application. Phys. Earth Planet. Inter. 183, 151–163.

Journaux, B., Brown, J.M., Pakhomova, A., Collings, I.E., Petitgirard, S., Espinoza, P., Boffa Ballaran, T., Vance, S.D., Ott, J., Cova, F., et al., 2020. Holistic approach for studying planetary hydrospheres: Gibbs representation of ices thermodynamics, elasticity, and the water phase diagram to 2,300 MPa. J. Geophy. Res. Planet. 125 (1), e2019JE006176.

Kamata, S., Nimmo, F., Sekine, Y., Kuramoto, K., Noguchi, N., Kimura, J., Tani, A., 2019. Pluto's ocean is capped and insulated by gas hydrates. Nat. Geosci. 12 (6), 407–410.

Kaula, W.M., 1964. Tidal dissipation by solid friction and the resulting orbital evolution. Rev. Geophys. 2 (4), 661–685.

Khan, A., Ceylan, S., van Driel, M., Giardini, D., Lognonné, P., Samuel, H., Schmerr, N.C., Stähler, S.C., Duran, A.C., Huang, Q., Kim, D., Broquet, A., Charalambous, C., Clinton, J.F., Davis, P.M., Drilleau, M., Karakostas, F., Lekic, V., McLennan, S.M., Maguire, R.R., Michaut, C., Panning, M.P., Pike, W.T., Pinot, B., Plasman, M., Scholz, J.-R., Widmer-Schnidrig, R., Spohn, T., Smrekar, S.E., Banerdt, W.B., 2021. Upper mantle structure of Mars from InSight seismic data. Science 373 (6553), 434–438.

Khan, A., Liebske, C., Rozel, A., Rivoldini, A., Nimmo, F., Connolly, J.A.D., Plesa, A.-C., Giardini, D., 2018. A geophysical perspective on the bulk composition of Mars. J. Geophys. Res. Planets 123 (2), 575–611.

Knapmeyer-Endrun, B., Panning, M.P., Bissig, F., Joshi, R., Khan, A., Kim, D., Lekić, V., Tauzin, B., Tharimena, S., Plasman, M., Compaire, N., Garcia, R.F., Margerin, L., Schimmel, M., Stutzmann, E., Schmerr, N., Bozdağ, E., Plesa, A.-C., Wieczorek, M.A., Broquet, A., Antonangeli, D., McLennan, S.M., Samuel, H., Michaut, C., Pan, L., Smrekar, S.E., Johnson, C.L., Brinkman, N., Mittelholz, A., Rivoldini, A., Davis, P.M., Lognonné, P., Pinot, B., Scholz, J.-R., Stähler, S., Knapmeyer, M., van Driel, M., Giardini, D., Banerdt, W.B., 2021. Thickness and structure of the martian crust from InSight seismic data. Science 373 (6553), 438–443.

Kovach, R.L., Chyba, C.F., 2001. Seismic detectability of a subsurface ocean on Europa. Icarus 150 (2), 279–287.

Lainey, V., Casajus, L.G., Fuller, J., Zannoni, M., Tortora, P., Cooper, N., Murray, C., Modenini, D., Park, R.S., Robert, V., et al., 2020. Resonance locking in giant planets indicated by the rapid orbital expansion of Titan. Nat. Astron. 4 (11), 1053–1058.

Lau, H.C., Faul, U.H., 2019. Elastostatic from seismic to tidal timescales: Theory and observations. Earth Planet. Sci. Lett. 508, 18–29.

Lee, S., Zanolin, M., Thode, A.M., Pappalardo, R.T., Makris, N.C., 2003. Probing Europa's interior with natural sound sources. Icarus 165 (1), 144–167.

Lodders, K., 2003. Solar system abundances and condensation temperatures of the elements. Astrophys. J. 591 (2), 1220.

Maguire, R.R., Schmerr, N.C., Lekić, V., Hurford, T.A., Dai, L., Rhoden, A.R., 2021. Constraining Europa's ice shell thickness with fundamental mode surface wave dispersion. Icarus 369, 114617.

Makarov, V.V., Efraimsky, M., 2013. No pseudosynchronous rotation for terrestrial planets and moons. Astrophys. J. 764 (1), 27.

Malhotra, R., Williams, J.G., 1997. Pluto's heliocentric orbit. Pluto Charon 127–157.

McCarthy, C., Cooper, R.F., 2016. Tidal dissipation in creeping ice and the thermal evolution of Europa. Earth Planet. Sci. Lett. 443, 185–194.

McCord, T.B., Sotin, C., 2005. Ceres: Evolution and current state. J. Geophy. Res. Planet. 110 (E5).

McKinnon, W.B., 2006. On convection in ice I shells of outer solar system bodies, with detailed application to callisto. Icarus 183 (2), 435–450.

McKinnon, W.B., Prialnik, D., Stern, S.A., Coradini, A., 2008. Structure and evolution of Kuiper belt objects and dwarf planets. Solar Syst. Beyond Neptune 1, 213–241.

McKinnon, W.B., Simonelli, D.P., Schubert, G., 1997. Composition, internal structure, and thermal evolution of Pluto and Charon. In: Stern, S.A., Tholen, D.J. (Eds.), Pluto And Charon. p. 295.

McKinnon, W.B., Stern, S., Weaver, H., Nimmo, F., Bierson, C., Grundy, W., Cook, J., Cruikshank, D., Parker, A., Moore, J., et al., 2017. Origin of the Pluto-Charon system: Constraints from the new horizons flyby. Icarus 287, 2–11.

Miller, N., Fortney, J., Jackson, B., 2009. Inflating and deflating hot Jupiters: Coupled tidal and thermal evolution of known transiting planets. Astrophys. J. 702 (2), 1413.

Moore, J.M., McKinnon, W.B., Spencer, J.R., Howard, A.D., Schenk, P.M., Beyer, R.A., Nimmo, F., Singer, K.N., Umurhan, O.M., White, O.L., et al., 2016. The geology of Pluto and Charon through the eyes of new horizons. Science 351 (6279), 1284–1293.

Moresi, L.N., Solomatov, V.S., 1995. Numerical investigation of 2D convection with extremely large viscosity variations. Phys. Fluids 7, 2154–2162.

Mottl, M.J., Glazer, B.T., Kaiser, R.I., Meech, K.J., 2007. Water and astrobiology. Geochemistry 67 (4), 253–282.

Neveu, M., Desch, S., Shock, E., Glein, C., 2015. Prerequisites for explosive cryovolcanism on dwarf planet-class Kuiper belt objects. Icarus 246, 48–64.

Neveu, M., Rhoden, A.R., 2019. Evolution of Saturn's mid-sized moons. Nat. Astron. 3 (6), 543–552.

Nimmo, F., Hamilton, D., McKinnon, W., Schenk, P., Binzel, R.P., Bierson, C., Beyer, R., Moore, J., Stern, S., Weaver, H., et al., 2016. Reorientation of sputnik planitia implies a subsurface ocean on Pluto. Nature 540 (7631), 94–96.

Nimmo, F., Pappalardo, R., 2016. Ocean worlds in the outer solar system. J. Geophy. Res. Planet. 121 (8), 1378–1399.

Nimmo, F., Umurhan, O., Lisse, C.M., Bierson, C.J., Lauer, T.R., Buie, M.W., Throop, H.B., Kammer, J.A., Roberts, J.H., McKinnon, W.B., et al., 2017. Mean radius and shape of Pluto and Charon from new horizons images. Icarus 287, 12–29.

Olkin, C.B., Ennico, K., Spencer, J., 2017. The Pluto system after the new horizons flyby. Nat. Astron. 1 (10), 663–670.

Olkin, C., Wasserman, L., Franz, O., 2003. The mass ratio of Charon to Pluto from hubble space telescope astrometry with the fine guidance sensors. Icarus 164 (1), 254–259.
Panning, M., Lekic, V., Manga, M., Cammarano, F., Romanowicz, B., 2006. Long-period seismology on Europa: 2. Predicted seismic response. J. Geophy. Res. Planet. 111 (E12).

Panning, M.P., Stähler, S.C., Huang, H.-H., Vance, S.D., Kedar, S., Tsai, V.C., Pike, W.T., Lorenz, R.D., 2018. Expected seismicity and the seismic noise environment of Europa. J. Geophy. Res. Planet. 123 (1), 163–179.

Pappalardo, R., Vance, S., Bagenal, F., Bills, B., Blaney, D., Blankenship, D., Brinckerhoff, W., Connerney, J., Hand, K., Hoehler, T.M., et al., 2013. Science potential from a Europa lander. Astrobiology 13 (8), 740–773.

Pathoff, D.A., Kattenhorn, S.A., 2011. A fracture history on enceladus provides evidence for a global ocean. Geophys. Res. Lett. 38 (18).

Renaud, J.P., Henning, W.G., 2018. Increased tidal dissipation using advanced rheological models: Implications for io and tidally active exoplanets. Astron. J. 857 (2), 98.

Renaud, J.P., Henning, W.G., Saxena, P., Neveu, M., Bagheri, A., Mandell, A., Hurford, T., 2021. Tidal dissipation in dual-body, highly eccentric, and non-synchronously rotating systems: Applications to Pluto–Charon and the exoplanet TRAPPIST-1e. Planet. Sci. J. 2 (1), 4.

Rhoden, A.R., Henning, W., Hurford, T.A., Hamilton, D.P., 2015. The interior and orbital evolution of Charon as preserved in its geologic record. Icarus 246, 11–20.

Rhoden, A.R., Miltitzer, B., Huff, E.M., Hurford, T.A., Manga, M., Richards, M.A., 2010. Constraints on Europa's rotational dynamics from modeling of tidally-driven fractures. Icarus 210 (2), 770–784.

Rhoden, A.R., Skjetne, H.L., Henning, W.G., Hurford, T.A., Walsh, K.J., Stern, S., Olkin, C., Spencer, J., Weaver, H., Young, L., et al., 2020. Charon: A brief history of tides. J. Geophy. Res. Planet. 125 (7), e2020JE006449.

Roberts, J.H., Nimmo, F., 2008. Tidal heating and the long-term stability of a subsurface ocean on enceladus. Icarus 194, 675–689.

Robuchon, G., Nimmo, F., 2011. Thermal evolution of pluto and implications for surface tectonics and a subsurface ocean. Icarus 216 (2), 426–439.

Rovira-Navarro, M., Rieutord, M., Gerkema, T., Maas, L.R., van der Wal, W., Vermeersen, B., 2019. Do tidally-generated inertial waves heat the subsurface oceans of europa and enceladus? Icarus 321, 126–140.

Samuel, H., Lognonné, P., Panning, M., Lainey, V., 2019. The rheology and thermal history of Mars revealed by the orbital evolution of phobos. Nature 569 (7757), 523.

Saxena, P., Renaud, J.P., Henning, W.G., Jutzi, M., Hurford, T., 2018. Relevance of tidal heating on large TNOs. Icarus 302, 245–260.

Schubert, G., Hussmann, H., Lainey, V., Matson, D., McKinnon, W., Sohl, F., Sotin, C., Tobie, G., Turrini, D., Van Hoolst, T., 2010. Evolution of icy satellites. Space Sci. Rev. 153 (1), 447–484.

Schubert, G., Spohn, T., Reynolds, R.T., 1986. Thermal histories, compositions and internal structures of the moons of the solar system. In: Burns, J.A., Matthews, M.S. (Eds.), Satellites. pp. 224–292.

Sekine, Y., Genda, H., Kamata, S., Funatsu, T., 2017. The Charon-forming giant impact as a source of Pluto's dark equatorial regions. Nat. Astron. 1 (2), 1–6.

Shoji, D., Hussmann, H., Sohl, F., Kurita, K., 2014. Non-steady state tidal heating of enceladus. Icarus 235, 75–85.

Shoji, D., Kurita, K., 2014. Thermal-orbital coupled tidal heating and habitability of martian-sized extrasolar planets around M stars. Astrophys. J. 789 (1), 3.

Spencer, J., Beyer, R.A., Robbins, S.J., Singer, K.N., Nimmo, F., 2021. The geology and geophysics of Charon. Pluto Syst. After New Horiz. 395.

Spencer, J.R., Grundy, W.M., Nimmo, F., Young, L.A., 2020. Chapter 12 - the Pluto system after new horizons. In: Prialnik, D., Barucci, M.A., Young, L.A. (Eds.), The Trans-Neptunian Solar System. Elsevier, pp. 271–288.

Stähler, S.C., Khan, A., Banerdt, W.B., Lognonné, P., Giardini, D., Ceylan, S., Dril-leau, M., Duran, A.C., Garcia, R.F., Huang, Q., Kim, D., Lekic, V., Samuel, H., Schimmel, M., Schmerr, N., Sollberger, D., Stutzmann, E., Xu, Z., Antonangeli, D., Charalambous, C., Davis, P.M., Irving, J.C.E., Kawamura, T., Knapmeyer, M., Maguire, R., Marusiak, A.G., Panning, M.P., Perrin, C., Plesa, A.-C., Rivoldini, A., Schmelzbach, C., Zenhäusern, G., Beucler, E., Clinton, J., Dahmen, N., van Driel, M., Gudkova, T., Horleston, A., Pike, W.T., Plasman, M., Smrekar, S.E., 2021. Seismic detection of the martian core. Science 373 (6553), 443–448.

Stähler, S.C., Panning, M.P., Vance, S.D., Lorenz, R.D., van Driel, M., Nissen-Meyer, T., Kedar, S., 2018. Seismic wave propagation in icy ocean worlds. J. Geophy. Res. Planet. 123 (1), 206–232.

Stern, S.A., Bagenal, F., Ennico, K., Gladstone, G., Grundy, W., McKinnon, W., Moore, J., Olkin, C., Spencer, J., Weaver, H., et al., 2015. The Pluto system: Initial results from its exploration by new horizons. Science 350 (6258).

Stern, S.A., Grundy, W.M., McKinnon, W.B., Weaver, H.A., Young, L.A., 2018. The Pluto system after new horizons. Annu. Rev. Astron. Astrophys. 56, 357–392.

Stern, S., Weaver, H., Steffl, A., Mutchler, M., Merline, W., Buie, M., Young, E., Young, L., Spencer, J., 2006. A giant impact origin for pluto's small moons and satellite multiplicity in the Kuiper belt. Nature 439 (7079), 946–948.

Sundberg, M., Cooper, R., 2010. A composite viscoelastic model for incorporating grain boundary sliding and transient diffusion creep; correlating creep and attenuation responses for materials with a fine grain size. Phil. Mag. 90 (20), 2817–2840.

Tobie, G., Grasset, O., Dumoulin, C., Mocquet, A., 2019. Tidal response of rocky and ice-rich exoplanets. Astron. Astrophys. 630, A70.

Travis, B., Olson, P., 1994. Convection with internal sources and turbulence in the Earth's mantle. Geophys. J. Int. 118, 1–19.

Tyler, R.H., 2008. Strong ocean tidal flow and heating on moons of the outer planets. Nature 456 (7223), 770–772.

Tyler, R.H., 2020. Heating of enceladus due to the dissipation of ocean tides. Icarus 348, 113821.

Tyler, R.H., Henning, W.G., Hamilton, C.W., 2015. Tidal heating in a magma ocean within Jupiter's moon Io. Astrophys. J. Suppl. Ser. 218 (2), 22.

Vance, S., Harnmeijer, J., Kimura, J., Hussmann, H., DeMartin, B., Brown, J.M., 2007. Hydrothermal systems in small ocean planets. Astrobiology 7 (6), 987–1005.

Vance, S.D., Kedar, S., Panning, M.P., Stahler, S.C., Bills, B.G., Lorenz, R.D., Huang, H.-H., Pike, W.T., Castillo, J.C., Lognonne, P., et al., 2018. Vital signs: seismology of icy ocean worlds. Astrobiology 18 (1), 37–53.

Vilella, K., Choblet, G., Tsao, W.-E., Deschamps, F., 2020. Tidally heated convection and the occurrence of melting in icy satellites: Application to Europa. J. Geophy. Res. Planet. 125 (3), e2019JE006248.

Waite Jr., J.H., Lewis, W., Magee, B., Lunine, J., McKinnon, W., Glein, C., Mousis, O., Young, D., Brockwell, T., Westlake, J., et al., 2009. Liquid water on enceladus from observations of ammonia and 40 Ar in the plume. Nature 460 (7254), 487–490.

Williams, J.G., Boggs, D.H., 2015. Tides on the Moon: Theory and determination of dissipation. J. Geophy. Res. Planet. 120 (4), 689–724.

Zannoni, M., Hemingway, D., Casajus, L.G., Tortora, P., 2020. The gravity field and interior structure of dione. Icarus 345, 113713.