# The Tidal-Thermal Evolution of the Pluto-Charon System

Amirhossein Bagheri1 Amir Khan1,2 Frederic Deschamps3 Henri Samuel4 Mikhail
Kruglyakov5,6 Domenico Giardini1 1Institute of Geophysics, ETH Zürich, Zürich,
Switzerland  2Physik-Institut, University of Zürich, Zürich, Switzerland  3
Institute of Earth Sciences, Academia Sinica, 128 Academia Road, Sector 2,
Nangang, Taipei 11529, Taiwan  4Institut de Physique du Globe de Paris, CNRS,
Université de Paris, Paris, France  5University of Otago, Dunedin, New Zealand
6Geoelectromagnetic Research Center, Institute of Physics of the Earth, Moscow

###### 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, including initial orbital and thermal parameters. 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 close to 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.

###### keywords:

Pluto Charon Tidal Evolution Thermal Evolution Subsurface Ocean

##  1 Introduction

The Kuiper belt harbours 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 further
away until it reached its current synchronous state (Dobrovolskis et al.,,
1997).

Generally speaking, tidal evolution drives planetary systems towards
equilibrium states by damping their orbital eccentricity, and forcing the spin
rates towards 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, influences 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 thermal 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 of the objects (distance, eccentricity, and
spin and orbit rates), these all need to be studied within a single framework.

In the context of tidal dissipation and the associated generation of heat, the
New Horizons mission found indications that Pluto and Charon may harbour
subsurface oceans beneath their ice-covered surfaces (Nimmo et al.,, 2017;
Olkin et al.,, 2017). The possibility that both bodies are able to sustain a
liquid ocean beneath an icy cover further makes them 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 eccentric and non-synchronously rotating system (Bagheri et al.,,
2021; Renaud et al.,, 2021). Specifically, tidal models that truncate
eccentricity functions to $e^{2}$ (where $e$ is eccentricity) on eccentric
orbits ($>$0.1) impart changes in spin rate evolution, spin-orbit resonances,
and errors in heating rates that typically increase significantly for very
high eccentricity ($>$0.5), which is not observed when higher-order terms are
included. Moreover, while the aforementioned studies found the evolution to be
fast ($\leq$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 combining the
comprehensive tidal model of Bagheri et al., (2021) that incorporates a proper
viscoelastic description of the tidal response of icy worlds with the
parameterized thermal convection models of Deschamps and Vilella, (2021) (icy
crust) that takes into account the effect of impurities in the liquid ocean
and Samuel et al., (2019) (silicate core). 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.

##  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 differentiated (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 differentiated (Stern et al.,, 2015, 2018; Spencer et al.,, 2021).
The observations include surface spectra dominated by ices and lack of
compression in the surface geological record that would otherwise be expected,
were the interiors undifferentiated as a consequence of the formation of
deeper and denser high-pressure ice 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 expansion that occurred when an early ocean froze above a
differentiated interior (Stern et al.,, 2015; Moore et al.,, 2016; Beyer et
al.,, 2017).

On account of a substantial rock fraction combined with a relatively large
size, it has been estimated that 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 that an ocean would never develop (Robuchon 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). In line herewith,
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).

Parameter | Symbol | Pluto | Charon
---|---|---|---
Radius (m) | $R$ | 1.1883$\times$ 106 | 0.606$\times$ 106
Mass (kg) | $M$ | 1.328$\times$ 1022 | 1.603$\times$ 1021
Density (kg/m3) | $\rho$ | 1854 | 1701
Semi-Major Axis (m) | $a$ | – | 19.596$\times$ 106
Spin period (day) | $\theta$ | 6.387 | 6.387
Eccentricity | $e$ | – |  $\sim$0
Inclination | $i$ | – |  $\sim$0
Table 1: Properties of the Pluto-Charon system (Nimmo et al.,, 2017; Stern et
al.,, 2015)

Plausible present-day interior structure models for Pluto and Charon are shown
as cross-sections in Figure 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 be
part of the silicate part hereafter, because of 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 3000–3600 kg/m3, as consequence of which the radius of
the silicate part varies in the range of 820–890 km for Pluto and 395–430 km
for Charon, to satisfy the mean density of the bodies. In turn, this implies
an 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.

Parameter/unit | Symbol | Value
---|---|---
Ice |  |
Shear modulus (GPa) | $\mu_{ice}$ | 4.8
Bulk modulus (GPa) | $\kappa_{ice}$ | 10.5
Density (kg/$\rm m^{3}$) | $\rho_{ice}$ | 920
Water |  |
Shear modulus (GPa) | $\mu_{w}$ | 0
Bulk modulus (GPa) | $\kappa_{w}$ | 2.2
Density ( kg/$\rm m^{3}$) | $\rho_{w}$ | 1000
Silicate |  |
Shear modulus (GPa) | $\mu_{sil}$ | 60
Bulk modulus (GPa) | $\kappa_{sil}$ | 120
Density (kg/$\rm m^{3}$) | $\rho_{sil}$ | 3000–3600
Table 2: Properties of different layers for the reference model Figure 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.

###  2.2 Orbital properties

The observed orbital properties, including semi-major axis ($a$), eccentricity
($e$), inclination ($i$), and spin period ($\theta$) of Pluto and Charon are
given 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 active
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 and to the surface, either through radiation,
conduction or convection. The main heat sources in a binary system are:

  1. 1.

the impact heating associated with accretion during planet formation,

  2. 2.

the gravitational energy released during planetary differentiation,

  3. 3.

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

  4. 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).

As illustrated in Figure 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.

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 large (e.g., Durham
et al.,, 2010), convection in the outer ice shell may occur in the so-called
stagnant-lid regime, for which the amount of heat that can be transported to
the surface is reduced 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). 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}\bigg{[}\rho_{w}C_{w}\bigg{(}\frac{\partial
T_{bot}}{\partial r}-\frac{\partial T_{ad}}{\partial
r}\bigg{)}\frac{r_{bot}^{3}-r_{c}^{3}}{3}-\rho_{I}L_{I}r_{bot}^{2}\bigg{]}=r_{bot}^{2}\phi_{bot}-r_{c}^{2}\phi_{c},$
(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})\bigg{[}1-\frac{\alpha_{w}}{\rho_{w}C_{w}}\rho_{I}g(r-r_{bot})\bigg{]},$
(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: $\rm{}^{235}U$, $\rm{}^{238}U$, 232Th, and 40K. 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}\bigg{(}r^{2}k_{c}\frac{\partial
T_{c}}{\partial r}\bigg{)}+H(t),$ (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 intial
conditions are the initial thermal profiles for either body (shown in Figure
3). Values for all parameters are listed in Table 3.

Parameter (Unit) | Symbol | Value/expression
---|---|---
Ice Ih properties |  |
Surface temperature $\rm(K)$ | $T_{surf}$ | 40
Thermal expansion $\rm(1/K)$ | $\alpha_{I}$ | 1.56$\times$ 10-4
Thermal conductivity $\rm(W/m/K)$ | $k$ | 2.6
Heat capacity $\rm(J/kg)$ | $C_{p}$ | 7.037$T$ \+ 185
Thermal diffusivity $\rm(m^{2}/s)$ | $\kappa_{I}$ | $k/{(\rho_{I}C_{p})}$
Latent heat of fusion $\rm(kJ/kg)$ | $L_{I}$ | 284
Reference bulk viscosity $\rm(Pa.s)$ | $\eta_{ref}$ | 1014
Activation energy $\rm(kJ/mol)$ | $E^{*}$ | 60
Initial Thickness (Pluto) $\rm(km)$ | $D_{init}^{ice}$ | 298
Initial Thickness (Charon) $\rm(km)$ | $D_{init}^{ice}$ | 171
Liquid water properties |  |
Initial ammonia content $\rm(wt\%)$ | $A_{init}$ | 3
Thermal expansion $\rm(1/K)$ | $\alpha_{w}$ | 3$\times$10-4
Heat capacity $\rm(J/kg)$ | $C_{w}$ | 4180
Initial Thickness (Pluto) $\rm(km)$ | $D_{init}^{w}$ | 15
Initial Thickness (Charon) $\rm(km)$ | $D_{init}^{w}$ | 15
Silicate core properties |  |
Thermal conductivity $\rm(W/m/K)$ | $k_{c}$ | 2.4
Heat capacity $\rm(J/kg/K)$ | $C_{c}$ | 1100
Core radius (Pluto) $(\rm km)$ | $R_{c}$ | 875
Core radius (Charon) $(\rm km)$ | $R_{c}$ | 420
Density $(\rm kg/m^{3})$ | $\rho_{c}$ | 3100
Table 3: Thermal properties of different layers of the system

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 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}$ ¿ $\phi_{cond}$.

To relate the ice shell Rayleigh number, $Ra$, which measures the vigour 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 TD^{3}}{\eta_{m}\kappa_{I}},$ (4)
$\gamma=\frac{E^{*}\Delta T}{R_{gas}{T_{m}}^{2}},$ (5)
$\phi_{bot}=1.46\frac{k_{I}\Delta T}{Df^{2}}\frac{Ra^{0.27}}{\gamma^{1.21}},$
(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}HD^{2}}{k_{I}\Delta
T}\frac{\Delta T}{Ra^{0.25}},$ (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\bigg{[}\frac{E^{*}}{R_{gas}T_{ref}}\bigg{(}\frac{T_{ref}}{T_{m}}-1\bigg{)}\bigg{]},$
(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_{H2O,bot}$. 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 theses 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 behaviour.
Associated herewith is an increase in bulk viscosity, which reduces the vigour
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 CH4, N2, and CO ices (Grundy et al.,, 2016). In
addition, H2O and NH3 ices have also been detected that could have originated
in the interior Dalle Ore et al., (2019). Here, we consider ammonia ($\rm
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 $\rm 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 $\rm H_{2}O-NH_{3}$ system, is that only
water crystallises until the eutectic composition is reached (32.2 vol% NH3).
Impurities (NH3) remain in the ocean, whereby their concentrations increase
and, as a result, their impact on 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 $\sim$270 K for pure water and $\sim$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_{l}({\bf{r}})=\left(\frac{R}{r}\right)^{l+1}k_{l}\Phi_{l}({\bf{R}},{\bf{r}}^{*}),$
(9)

where $R$ is the radius of the planet, ${\bf{R}}$ is a point on the planet’s
surface, ${\bf{r}}$ is a point exterior to ${\bf{R}}$, while ${\bf{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_{\mathrm{l}}$
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 centre towards 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_{l}({\bf{r}},t)=\left(\frac{R}{r}\right)^{l+1}\hat{\mathrm{k}}_{l}\Phi_{l}({\bf{R}},{\bf{r}}^{*}),$
(10)

where $\hat{\mathrm{k}}_{l}$ is the (linear) Love operator that maps the
entire history of the perturbation ($\Phi_{l}(t^{\prime})$ over
$t^{\prime}\leq 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}({\bf{r}},\omega^{\mathrm{nm}}_{\mathrm{pq}})=\left(\frac{R}{r}\right)^{l+1}\bar{\mathrm{k}}_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})\bar{\Phi}_{l}({\bf{R}},{\bf{r}}^{*},\omega^{\mathrm{nm}}_{\mathrm{pq}}),$
(11)

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

$\bar{\mathrm{k}}_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})=\mathrm{Re}[\bar{\mathrm{k}}_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})]+i\,\mathrm{Im}[\bar{\mathrm{k}}_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})]=|\bar{\mathrm{k}}_{l}|~{}e^{-\epsilon_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})}.$
(12)

In Eq. 11, $\bar{\psi}_{l}({\bf{r}},\omega^{\mathrm{nm}}_{\mathrm{pq}})$ is
lagging behind
$\bar{\Phi}_{l}({\bf{R}},{\bf{r}}^{*},\omega^{\mathrm{nm}}_{\mathrm{pq}})$ by
the phase angle $\epsilon_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})$, which by
convention is the negative argument of the complex Love number
$\bar{\mathrm{k}}_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})$. Replacing
${\bf{r}}$ by the position of the perturbing body ${\bf{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{nm}}_{\mathrm{pq}}$ 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{nm}}_{\mathrm{pq}})}=\sin|\epsilon_{l}(\omega^{\mathrm{nm}}_{\mathrm{pq}})|.$
(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 more “flexibility” of Pluto
relative to 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 the 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, G. 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
600$\rm{}^{\circ}C$, viscoelastic creep does not occur (Jackson and Faul,,
2010). In the case of Pluto and Charon, temperatures of the silicate core do
not exceed 500$\rm{}^{\circ}C$ 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
further discussed 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 employ the tidal evolution model of Bagheri et al., (2021), which
extends the Darwin-Kaula 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 higher-order eccentricity functions and harmonic
modes is essential. Here, we further extended our tidal model to include the
case of non-synchronous rotation. Due to the spherical shape of both bodies,
librational tides are not important.

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^{\prime}$), the planet and satellite radii ($R$ and $R^{\prime}$), the
quality functions ($K_{l}=k_{l}/Q_{l}$ and
$K_{l}^{\,\prime}=k_{l}^{\,\prime}/Q_{l}^{\,\prime}$), and spin rates
($\dot{\theta}$, and $\dot{\theta}^{\,\prime}$):

$\frac{da}{dt}=\frac{G(M+M^{\prime})}{n~{}a^{2}}\frac{M^{\prime}}{M}\bigg{(}\frac{R}{a}\bigg{)}^{5}\bigg{[}\underbrace{{\cal{F}}(K_{l},\,\dot{\theta},\,n,\,e)}_{\text{tides
on primary
body}}+\underbrace{{\cal{F}}(K^{\prime}_{l},\,\dot{\theta}^{\prime},n,\,e)}_{\text{tides
on satellite body}}\bigg{]},$ (14)
$\frac{de}{dt}=\frac{G(M+M^{\prime})}{n~{}a^{3}}\frac{M^{\prime}}{M}\bigg{(}\frac{R}{a}\bigg{)}^{5}\bigg{[}\underbrace{{\cal{L}}(K_{l},\,\dot{\theta},\,n,\,e)}_{\text{tides
on
primary}}+\underbrace{{\cal{L}}(K^{\prime}_{l},\,\dot{\theta}^{\prime},n,\,e)}_{\text{tides
on satellite}}\bigg{]}\;\;,$ (15)
$\bigg{(}\frac{d\dot{\theta}}{dt}\bigg{)}_{\text{primary}}=~{}\frac{GM^{\prime
2}}{a^{3}MR^{2}}\bigg{(}\frac{R}{a}\bigg{)}^{3}\bigg{[}\underbrace{{\cal{G}}(K_{l},\,\dot{\theta},\,n,\,e)}_{\text{tides
on primary}}\bigg{]},\;\;$ (16)
$\bigg{(}\frac{di}{dt}\bigg{)}_{\text{primary}}=~{}n\sin
i\frac{M^{\prime}}{M}\bigg{(}\frac{R}{a}\bigg{)}^{5}\bigg{[}\underbrace{{{\cal{I}}(K_{l},\dot{\theta},n,i,e)}}_{\text{tides
on primary}}\bigg{]},\;\ $ (17)
$\bigg{(}\frac{dE}{dt}\bigg{)}_{\text{primary}}=~{}\frac{GM^{\prime
2}}{a}\bigg{(}\frac{R}{a}\bigg{)}^{5}\bigg{[}\underbrace{{\cal{H}}(K_{l},\,\dot{\theta},\,n,\,e)}_{\text{tides
on primary}}\bigg{]}.\;\;$ (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^{\prime}$, $K_{l}^{\,\prime}$, and
$\dot{\theta}^{\prime}$, and interchanging $M$ and $M^{\prime}$. Detailed
expressions for the functions ${\cal{F}}$, ${\cal{L}}$, ${\cal{I}}$,
${\cal{G}}$, and ${\cal{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 employed are $e$$\sim$0.4,
$a$$\sim$0.65$a_{p}$, and $\sim$5 times the initial mean motion, corresponding
to rotation period around 10 hr, respectively, where $a_{p}$ is the present-
day observed semi-major axis. These values 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. The choice of initial parameters
will be discussed further in section 4.3.

###  4.1 Tidal evolution

Figure 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. Figure 2a shows that the tidal evolution of the binary system
occurs very rapidly with the stable 1:1 spin-orbit resonance state
($\dot{\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 ($\dot{\theta}/n=1.5$) after $\sim$104 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.

Figure 2b shows the expansion of the semi-major axis. As illustrated in the
plot, the orbit does not expand monotonically because of the aforementioned
resonances. In fact, 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. Once the satellite falls out of this resonance and has
reached the 1:1 spin-orbit resonance, Pluto is propelled toward
synchronisation with the orbital motion $n$. 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.

Related to the spin rate and semi-major axis behaviour is the change in
eccentricity. Figure 2c shows that the eccentricity of the orbit increases for
a very short time in the beginning before Charon’s fall into the 3:2 spin-
orbit resonance, at which point, the eccentricity starts to dampen smoothly to
0, dictating the fall of Charon into the tidally-locked stated. A
qualitatively similar behaviour in the orbital evolution was observed by
Renaud et al., (2021), but because of a different initial $\dot{\theta}/n$ and
differences in planetary structure (physical properties and attenuation), the
exact timing of the 3:2 and 1:1 resonance captures vary slightly in comparison
to our results.

Finally, Figure 2d displays the heat generated in the two bodies due to both
tidal and radiogenic heating. We can make two 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, but that the tidal heat
source rapidly decreases as the eccentricity of the orbit diminishes; 3)
radiogenic heating appears to govern the thermal evolution throughout most of
the history of both bodies; and 4) 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 can potentially manifest in the form of surface features
that survive until the present time, and can be used to constrain the past
history of the objects (e.g., Robuchon and Nimmo,, 2011; Rhoden et al.,,
2020).

Figure 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_{p}$ are the mean anomaly and the
present-day distance between Pluto and Charon listed in Table 1.

###  4.2 Evolution of thermal structure

The thermal evolution of Pluto and Charon for the nominal case is shown in
Figure 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 Figure 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.

Figure 3 illustrates the thermal evolution of Pluto (Figure 3a) and Charon
(Figure 3b) with time in terms of snapshots of 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 $\sim$1300 K (Pluto) and
$\sim$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 104 yr (Figure 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 (Figure 3c) on account
of its size and therefore limited radiogenic budget. As illustrated in Figure
3a (inset), convection on Pluto starts operating in the ice shell after
$\sim$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
effectively removed through conduction, and convection in the bottom part of
the ice layer never sets in (Bierson et al.,, 2018; Nimmo et al.,, 2017). In
this particular simulation, Charon’s subsurface ocean solidified completely
$\sim$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 H2O
crystallizes out of the ocean. As a consequence, the concentration of NH3 in
the ocean grows, reaching the eutectic at $\sim$3.5 Gyr, at which point we
observe a change in slope in the crystallization behaviour of the ocean, which
is due to the fact that both NH3 and H2O start to crystallize and the
concentration of NH3 in the remaining ocean no longer changes (on Pluto the
eutectic is never reached). An interesting consequence of this result is that
the bottom of the ice shell on Charon may be composed of a mix of H2O and NH3,
assuming that the latter compound was present in the initial ocean.

Figure 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 $\sim$110 km and $\sim$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 $\sim$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 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.

Figure 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.

###  4.3 Influence of model parameters

####  4.3.1 Orbital parameters

In the nominal case (Figure 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 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 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 to arrive at its
final state and the amount of tidally-dissipated heat.

As part of our simulations, we also considered cases where $a_{0}>a_{p}$,
i.e., initial planet-satellite separation greater than the 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$106 years and, relative to the contribution
from radioactive decay, tidal heating has little affect 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 $\sim$3.5 Gyr, has refrozen completely (Figure 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 (Figure 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 variations with temperature.
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)). Figure 4a displays the
evolution of ocean thickness on Pluto and Charon for several different values
of reference viscosity ($\eta_{ref}$): 1012, 1014, and 1016 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
$\sim$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 $\sim$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 $\sim$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. 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 $\sim$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/cm3) to ensure that
total mass of the body remains unchanged. All other relevant parameters are as
in Table 3. The results are shown in Figure 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 ($\leq$5 km) subsurface ocean on Charon at the present day.
Moreover, the changes in core size result in subsurface oceans that vary in
thickness by $\sim$30 km and $\sim$20 km on Pluto and Charon, respectively.
Again, convection on Charon is generally inhibited by the presence of a
smaller hydrosphere.

We further investigated the effect of initial ocean/ice layer thickness
(Figure 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 (Figure 4a) for the low reference viscosity
cases ($\eta_{ref}$=1012 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, Figure 4c
shows that subsurface oceans on Pluto and Charon reach thicknesses $\sim$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.

Figure 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 ($\sim$3.5 Gyr) is indicative of a change in the crystallization
behaviour of the ocean. See main text for more details.

##### 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
Figure 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
$\sim$3.5–4 Gyr. Between the various cases, ocean thicknesses on Pluto and
Charon range from $\sim$10–50 km and $\sim$40–120 km, respectively. Convection
only happens for NH3=1 and 3 wt%, respectively, whereas for NH3=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 NH3=1 wt%,
convection is only operative between $\sim$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 Pluto and Charon.

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 $\sim$3.5 Gyr), but because the ice layer
is relatively thin, convection is unable to operate effectively, as a
consequence of which 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. Observations of extensional and compressional surface
tectonics on Charon and Pluto, respectively, made by New Horizons supports
this. While we applied our thermal-tidal 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 $\sim$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). 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. 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, for example.

##  6 Acknowledgements

This work was supported by a grant from the Swiss National Science Foundation
(project 172508). We would like to acknowledge Joe Renaud and Michael
Efroimsky for informed discussions on tidal evolution and dissipation.

## References

  * Anderson and Suga, (1994) Anderson, O. and Suga, H. (1994).  Thermal conductivity of the ih and xi phases of ice.  Physical Review B, 50:6583–6588.
  * Arakawa et al., (2019) Arakawa, S., Hyodo, R., and Genda, H. (2019).  Early formation of moons around large trans-neptunian objects via giant impacts.  Nature Astronomy, 3(9):802–807.
  * Asphaug and Reufer, (2013) Asphaug, E. and Reufer, A. (2013).  Late origin of the saturn system.  Icarus, 223(1):544–565.
  * Bagheri et al., (2019) Bagheri, A., Khan, A., Al-Attar, D., Crawford, O., and Giardini, D. (2019).  Tidal response of Mars constrained from laboratory-based viscoelastic dissipation models and geophysical data.  J. Geophys. Res. planets.
  * Bagheri et al., (2021) Bagheri, A., Khan, A., Efroimsky, M., Kruglyakov, M., and Giardini, D. (2021).  Dynamical evidence for phobos and deimos as remnants of a disrupted common progenitor.  Nature Astronomy, pages 1–5.
  * Barnes et al., (2021) Barnes, J. W., Turtle, E. P., Trainer, M. G., Lorenz, R. D., MacKenzie, S. M., Brinckerhoff, 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., and Stähler, S. C. (2021).  Science goals and objectives for the dragonfly titan rotorcraft relocatable lander.  The Planetary Science Journal, 2(4):130.
  * Barr and Collins, (2015) Barr, A. C. and Collins, G. C. (2015).  Tectonic activity on pluto after the charon-forming impact.  Icarus, 246:146–155.
  * Beuthe, (2019) Beuthe, M. (2019).  Enceladus’s crust as a non-uniform thin shell: Ii tidal dissipation.  Icarus, 332:66–91.
  * Beyer et al., (2017) 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., and Zangari, A. M. (2017).  Charon tectonics.  Icarus, 287:161–174.  Special Issue: The Pluto System.
  * Beyer et al., (2019) 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 et al., (2018) Bierson, C., Nimmo, F., and McKinnon, W. (2018).  Implications of the observed pluto–charon density contrast.  Icarus, 309:207–219.
  * Bierson et al., (2020) Bierson, C. J., Nimmo, F., and Stern, S. A. (2020).  Evidence for a hot start and early ocean formation on pluto.  Nature Geoscience, 13(7):468–472.
  * Bierson and Steinbrügge, (2021) Bierson, C. J. and Steinbrügge, G. (2021).  Tidal heating did not dry out io and europa.  The Planetary Science Journal, 2(3):89.
  * Boué and Efroimsky, (2019) Boué, G. and Efroimsky, M. (2019).  Tidal evolution of the Keplerian elements.  Celest. Mech. Dyn. Astron., 131(7):30.
  * Brozović et al., (2015) Brozović, M., Showalter, M. R., Jacobson, R. A., and Buie, M. W. (2015).  The orbits and masses of satellites of pluto.  Icarus, 246:317–329.
  * Buie et al., (2012) Buie, M. W., Tholen, D. J., and Grundy, W. M. (2012).  The orbit of charon is circular.  The Astronomical Journal, 144(1):15.
  * Canup and Salmon, (2018) Canup, R. and 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.
  * Canup, (2005) Canup, R. M. (2005).  A giant impact origin of pluto-charon.  Science, 307(5709):546–550.
  * Canup, (2010) Canup, R. M. (2010).  On a giant impact origin of charon, nix, and hydra.  The Astronomical Journal, 141(2):35.
  * Canup et al., (2021) Canup, R. M., Kratter, K. M., and Neveu, M. (2021).  On the origin of the pluto system.  The Pluto System After New Horizons, page 475.
  * Caswell and Cooper, (2016) Caswell, T. and Cooper, R. (2016).  Convection-induced microstructure and tidal dissipation in polycrystalline ice; an experimental approach.  In LPSC, number 1903 in Lunar and Planetary Science Conference, page 1129.
  * Caswell et al., (2015) Caswell, T. E., Cooper, R. F., and Goldsby, D. L. (2015).  The constant-hardness creep compliance of polycrystalline ice.  Geophysical Research Letters, 42(15):6261–6268.
  * Cheng et al., (2014) Cheng, W., Lee, M. H., and Peale, S. (2014).  Complete tidal evolution of pluto–charon.  Icarus, 233:242–258.
  * Citron et al., (2015) Citron, R. I., Genda, H., and Ida, S. (2015).  Formation of Phobos and Deimos via a giant impact.  Icarus, 252:334–338.
  * Clark et al., (2014) Clark, R. N., Swayze, G. A., Carlson, R., Grundy, W., and Noll, K. (2014).  Spectroscopy from space.  Reviews in Mineralogy and Geochemistry, 78(1):399–446.
  * Dalle Ore et al., (2019) 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., Verbiscer, A., Parker, A. H., Singer, K. N., Umurhan, O. M., Weaver, H. A., Olkin, C. B., Young, L. A., and Ennico, K. (2019).  Detection of ammonia on pluto’s surface in a region of geologically recent tectonism.  Science Advances, 5(5).
  * Dalle Ore et al., (2018) Dalle Ore, C. M., Protopapa, S., Cook, J., Grundy, W., Cruikshank, D., Verbiscer, 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 and Jaupart, (1993) Davaille, A. and Jaupart, C. (1993).  Transient high-rayleigh-number thermal convection with large viscosity variations.  J. Fluid Mech., 253:141–166.
  * Desch and Neveu, (2017) Desch, S. and Neveu, M. (2017).  Differentiation and cryovolcanism on charon: A view before and after new horizons.  Icarus, 287:175–186.
  * Desch, (2015) Desch, S. J. (2015).  Density of charon formed from a disk generated by the impact of partially differentiated bodies.  Icarus, 246:37–47.
  * Desch et al., (2009) Desch, S. J., Cook, J. C., Doggett, T., and Porter, S. B. (2009).  Thermal evolution of kuiper belt objects, with implications for cryovolcanism.  Icarus, 202(2):694–714.
  * Deschamps, (2021) Deschamps, F. (2021).  Stagnant lid convection with temperature-dependent thermal conductivity and the thermal evolution of icy worlds.  Geophysical Journal International, 224(3):1870–1889.
  * Deschamps and Sotin, (2001) Deschamps, F. and Sotin, C. (2001).  Thermal convection in the outer shell of large icy satellites.  Journal of Geophysical Research: Planets, 106(E3):5107–5121.
  * Deschamps et al., (2010) Deschamps, F., Tackley, P. J., and Nakagawa, T. (2010).  Temperature and heat flux scalings for isoviscous thermal convection in spherical geometry.  Geophys. J. Int., 182:137–154.
  * Deschamps and Vilella, (2021) Deschamps, F. and Vilella, K. (2021).  Scaling laws for mixed-heated stagnant lid convection and application to europa.  JGR Planets, in press.
  * Dobos et al., (2019) Dobos, V., Barr, A. C., and Kiss, L. L. (2019).  Tidal heating and the habitability of the trappist-1 exoplanets.  Astronomy & Astrophysics, 624:A2.
  * Dobrovolskis et al., (1997) Dobrovolskis, A., Peale, S., and Harris, A. (1997).  Dynamics of the pluto-charon binary.  Pluto and Charon, page 159.
  * Driscoll and Barnes, (2015) Driscoll, P. E. and Barnes, R. (2015).  Tidal heating of earth-like exoplanets around m stars: thermal, magnetic, and orbital evolutions.  Astrobiology, 15(9):739–760.
  * Durham et al., (2010) Durham, W. B., Prieto-Ballesteros, O., Goldsby, D. L., and Kargel, J. S. (2010).  Rheological and thermal properties of icy minerals.  Space Sci. Rev., 153:273–298.
  * Efroimsky, (2012) Efroimsky, M. (2012).  Tidal dissipation compared to seismic dissipation: In small bodies, Earths, and super-Earths.  Astrophys. J., 746:150.
  * Efroimsky, (2015) Efroimsky, M. (2015).  Tidal evolution of asteroidal binaries. ruled by viscosity. ignorant of rigidity.  Astron. J., 150(4):98.
  * Efroimsky and Makarov, (2013) Efroimsky, M. and Makarov, V. V. (2013).  Tidal friction and tidal lagging. applicability limitations of a popular formula for the tidal torque.  Astron. J., 764:26.
  * Efroimsky and Makarov, (2014) Efroimsky, M. and Makarov, V. V. (2014).  Tidal dissipation in a homogeneous spherical body. i. methods.  Astron. J., 795(1):6.
  * Egbert and Ray, (2003) Egbert, G. D. and Ray, R. D. (2003).  Semi-diurnal and diurnal tidal dissipation from topex/poseidon altimetry.  Geophysical Research Letters, 30(17).
  * Ferraz-Mello et al., (2020) Ferraz-Mello, S., Beaugé, C., Folonier, H. A., and Gomes, G. O. (2020).  Tidal friction in satellites and planets. the new version of the creep tide theory.  The European Physical Journal Special Topics, 229:1441–1462.
  * Grasset and Pargamin, (2005) Grasset, O. and Pargamin, J. (2005).  The ammonia–water system at high pressures: Implications for the methane of titan.  Planetary and Space Science, 53(4):371–384.
  * Grasset and Sotin, (1996) Grasset, O. and Sotin, C. (1996).  The cooling rate of a liquid shell in titan’s interior.  Icarus, 123(1):101–112.
  * Grundy et al., (2016) 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 et al., (2016) Hammond, N. P., Barr, A. C., and Parmentier, E. M. (2016).  Recent tectonic activity on pluto driven by phase changes in the ice shell.  Geophysical Research Letters, 43(13):6775–6782.
  * Harada et al., (2014) Harada, Y., Goossens, S., Matsumoto, K., Yan, J., Ping, J., Noda, H., and 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 et al., (2011) Heller, R., Leconte, J., and Barnes, R. (2011).  Tidal obliquity evolution of potentially habitable planets.  Astron. Astrophys., 528:A27.
  * Hussmann et al., (2010) Hussmann, H., Choblet, G., Lainey, V., Matson, D. L., Sotin, C., Tobie, G., and Van Hoolst, T. (2010).  Implications of rotation, orbital states, energy sources, and heat transport for internal processes in icy satellites.  Space science reviews, 153(1):317–348.
  * Hussmann et al., (2006) Hussmann, H., Sohl, F., and 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 et al., (2007) Hussmann, H., Sotin, C., and Lunine, J. I. (2007).  Interiors and evolution of icy satellites, in planets and moons.  Treatise on Geophysics, 10:509–539.
  * Hussmann and Spohn, (2004) Hussmann, H. and Spohn, T. (2004).  Thermal-orbital evolution of io and europa.  Icarus, 171(2):391 – 410.
  * Ida et al., (1997) Ida, S., Canup, R. M., and Stewart, G. R. (1997).  Lunar accretion from an impact-generated disk.  Nature, 389(6649):353–357.
  * Jackson, (2005) 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, pages 95–119. Elsevier.
  * Jackson and Faul, (2010) Jackson, I. and 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.
  * Kamata et al., (2019) Kamata, S., Nimmo, F., Sekine, Y., Kuramoto, K., Noguchi, N., Kimura, J., and Tani, A. (2019).  Pluto’s ocean is capped and insulated by gas hydrates.  Nature Geoscience, 12(6):407–410.
  * Khan et al., (2021) 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., and Banerdt, W. B. (2021).  Upper mantle structure of mars from insight seismic data.  Science, 373(6553):434–438.
  * Khan et al., (2018) Khan, A., Liebske, C., Rozel, A., Rivoldini, A., Nimmo, F., Connolly, J. A. D., Plesa, A., and Giardini, D. (2018).  A geophysical perspective on the bulk composition of Mars.  J. Geophys. Res. Planets, 123(2):575–611.
  * Knapmeyer-Endrun et al., (2021) 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, É., 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., and Banerdt, W. B. (2021).  Thickness and structure of the martian crust from insight seismic data.  Science, 373(6553):438–443.
  * Kovach and Chyba, (2001) Kovach, R. L. and Chyba, C. F. (2001).  Seismic detectability of a subsurface ocean on europa.  Icarus, 150(2):279–287.
  * Lainey et al., (2020) 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.  Nature Astronomy, 4(11):1053–1058.
  * Lau and Faul, (2019) Lau, H. C. and Faul, U. H. (2019).  Anelasticity from seismic to tidal timescales: Theory and observations.  Earth and Planetary Science Letters, 508:18–29.
  * Lee et al., (2003) Lee, S., Zanolin, M., Thode, A. M., Pappalardo, R. T., and Makris, N. C. (2003).  Probing europa’s interior with natural sound sources.  Icarus, 165(1):144–167.
  * Lodders, (2003) Lodders, K. (2003).  Solar system abundances and condensation temperatures of the elements.  The Astrophysical Journal, 591(2):1220.
  * Maguire et al., (2021) Maguire, R. R., Schmerr, N. C., Lekić, V., Hurford, T. A., Dai, L., and Rhoden, A. R. (2021).  Constraining europa’s ice shell thickness with fundamental mode surface wave dispersion.  Icarus, 369:114617.
  * Makarov and Efroimsky, (2013) Makarov, V. V. and Efroimsky, M. (2013).  No pseudosynchronous rotation for terrestrial planets and moons.  The Astrophysical Journal, 764(1):27.
  * Malhotra and Williams, (1997) Malhotra, R. and Williams, J. G. (1997).  Pluto’s heliocentric orbit.  Pluto and Charon, pages 127–157.
  * McCarthy and Cooper, (2016) McCarthy, C. and Cooper, R. F. (2016).  Tidal dissipation in creeping ice and the thermal evolution of europa.  Earth and Planetary Science Letters, 443:185–194.
  * McKinnon et al., (1997) McKinnon, W., Simonelli, D., and Schubert, G. (1997).  Composition, internal structure, and thermal evolution of pluto and charon.  Pluto and Charon, page 295.
  * McKinnon, (2006) 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 et al., (2008) McKinnon, W. B., Prialnik, D., Stern, S. A., and Coradini, A. (2008).  Structure and evolution of kuiper belt objects and dwarf planets.  The solar system beyond Neptune, 1:213–241.
  * McKinnon et al., (2017) 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 et al., (2009) Miller, N., Fortney, J., and Jackson, B. (2009).  Inflating and deflating hot jupiters: Coupled tidal and thermal evolution of known transiting planets.  The Astrophysical Journal, 702(2):1413.
  * Moore et al., (2016) 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 and Solomatov, (1995) Moresi, L. N. and Solomatov, V. S. (1995).  Numerical investigation of 2d convection with extremely large viscosity variations.  Phys. Fluids, 7:2154–2162.
  * Mottl et al., (2007) Mottl, M. J., Glazer, B. T., Kaiser, R. I., and Meech, K. J. (2007).  Water and astrobiology.  Geochemistry, 67(4):253–282.
  * Neveu et al., (2015) Neveu, M., Desch, S., Shock, E., and Glein, C. (2015).  Prerequisites for explosive cryovolcanism on dwarf planet-class kuiper belt objects.  Icarus, 246:48–64.
  * Neveu and Rhoden, (2019) Neveu, M. and Rhoden, A. R. (2019).  Evolution of saturn’s mid-sized moons.  Nature astronomy, 3(6):543–552.
  * Nimmo and Pappalardo, (2016) Nimmo, F. and Pappalardo, R. (2016).  Ocean worlds in the outer solar system.  Journal of Geophysical Research: Planets, 121(8):1378–1399.
  * Nimmo et al., (2017) 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 et al., (2003) Olkin, C., Wasserman, L., and 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.
  * Olkin et al., (2017) Olkin, C. B., Ennico, K., and Spencer, J. (2017).  The pluto system after the new horizons flyby.  Nature Astronomy, 1(10):663–670.
  * Panning et al., (2006) Panning, M., Lekic, V., Manga, M., Cammarano, F., and Romanowicz, B. (2006).  Long-period seismology on europa: 2. predicted seismic response.  Journal of Geophysical Research: Planets, 111(E12).
  * Panning et al., (2018) Panning, M. P., Stähler, S. C., Huang, H.-H., Vance, S. D., Kedar, S., Tsai, V. C., Pike, W. T., and Lorenz, R. D. (2018).  Expected seismicity and the seismic noise environment of europa.  Journal of Geophysical Research: Planets, 123(1):163–179.
  * Pappalardo et al., (2013) 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.
  * Renaud and Henning, (2018) Renaud, J. P. and Henning, W. G. (2018).  Increased tidal dissipation using advanced rheological models: Implications for io and tidally active exoplanets.  Astron. J., 857(2):98.
  * Renaud et al., (2021) Renaud, J. P., Henning, W. G., Saxena, P., Neveu, M., Bagheri, A., Mandell, A., and Hurford, T. (2021).  Tidal dissipation in dual-body, highly eccentric, and nonsynchronously rotating systems: Applications to pluto–charon and the exoplanet trappist-1e.  The Planetary Science Journal, 2(1):4.
  * Rhoden et al., (2015) Rhoden, A. R., Henning, W., Hurford, T. A., and Hamilton, D. P. (2015).  The interior and orbital evolution of charon as preserved in its geologic record.  Icarus, 246:11–20.
  * Rhoden et al., (2020) 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.  Journal of Geophysical Research: Planets, 125(7):e2020JE006449.
  * Roberts and Nimmo, (2008) Roberts, J. H. and Nimmo, F. (2008).  Tidal heating and the long-term stability of a subsurface ocean on Enceladus.  Icarus, 194:675–689.
  * Robuchon and Nimmo, (2011) Robuchon, G. and Nimmo, F. (2011).  Thermal evolution of pluto and implications for surface tectonics and a subsurface ocean.  Icarus, 216(2):426–439.
  * Samuel et al., (2019) Samuel, H., Lognonné, P., Panning, M., and Lainey, V. (2019).  The rheology and thermal history of Mars revealed by the orbital evolution of Phobos.  Nature, 569(7757):523.
  * Saxena et al., (2018) Saxena, P., Renaud, J. P., Henning, W. G., Jutzi, M., and Hurford, T. (2018).  Relevance of tidal heating on large tnos.  Icarus, 302:245–260.
  * Schubert et al., (2010) Schubert, G., Hussmann, H., Lainey, V., Matson, D., McKinnon, W., Sohl, F., Sotin, C., Tobie, G., Turrini, D., and Van Hoolst, T. (2010).  Evolution of icy satellites.  Space Science Reviews, 153(1):447–484.
  * Schubert et al., (1986) Schubert, G., Spohn, T., and Reynolds, R. T. (1986).  Thermal histories, compositions and internal structures of the moons of the solar system., pages 224–292.  University of Arizona Press.
  * Sekine et al., (2017) Sekine, Y., Genda, H., Kamata, S., and Funatsu, T. (2017).  The charon-forming giant impact as a source of pluto’s dark equatorial regions.  Nature Astronomy, 1(2):1–6.
  * Shoji et al., (2014) Shoji, D., Hussmann, H., Sohl, F., and Kurita, K. (2014).  Non-steady state tidal heating of enceladus.  Icarus, 235:75–85.
  * Shoji and Kurita, (2014) Shoji, D. and Kurita, K. (2014).  Thermal-orbital coupled tidal heating and habitability of martian-sized extrasolar planets around m stars.  The Astrophysical Journal, 789(1):3.
  * Spencer et al., (2021) Spencer, J., Beyer, R. A., Robbins, S. J., Singer, K. N., and Nimmo, F. (2021).  The geology and geophysics of charon.  The Pluto System After New Horizons, page 395.
  * Spencer et al., (2020) Spencer, J. R., Grundy, W. M., Nimmo, F., and Young, L. A. (2020).  Chapter 12 - the pluto system after new horizons.  In Prialnik, D., Barucci, M. A., and Young, L. A., editors, The Trans-Neptunian Solar System, pages 271–288. Elsevier.
  * Stähler et al., (2021) Stähler, S. C., Khan, A., Banerdt, W. B., Lognonné, P., Giardini, D., Ceylan, S., Drilleau, M., Duran, A. C., Garcia, R. F., Huang, Q., Kim, D., Lekic, V., Samuel, H., Schimmel, M., Schmerr, N., Sollberger, D., Stutzmann, É., 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, É., Clinton, J., Dahmen, N., van Driel, M., Gudkova, T., Horleston, A., Pike, W. T., Plasman, M., and Smrekar, S. E. (2021).  Seismic detection of the martian core.  Science, 373(6553):443–448.
  * Stähler et al., (2018) Stähler, S. C., Panning, M. P., Vance, S. D., Lorenz, R. D., van Driel, M., Nissen-Meyer, T., and Kedar, S. (2018).  Seismic wave propagation in icy ocean worlds.  Journal of Geophysical Research: Planets, 123(1):206–232.
  * Stern et al., (2006) Stern, S., Weaver, H., Steffl, A., Mutchler, M., Merline, W., Buie, M., Young, E., Young, L., and Spencer, J. (2006).  A giant impact origin for pluto’s small moons and satellite multiplicity in the kuiper belt.  Nature, 439(7079):946–948.
  * Stern et al., (2015) 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 et al., (2018) Stern, S. A., Grundy, W. M., McKinnon, W. B., Weaver, H. A., and Young, L. A. (2018).  The pluto system after new horizons.  Annual Review of Astronomy and Astrophysics, 56:357–392.
  * Sundberg and Cooper, (2010) Sundberg, M. and 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.  Philos. Mag., 90(20):2817–2840.
  * Tobie, G. et al., (2019) Tobie, G., Grasset, O., Dumoulin, C., and Mocquet, A. (2019).  Tidal response of rocky and ice-rich exoplanets.  A&A, 630:A70.
  * Travis and Olson, (1994) Travis, B. and Olson, P. (1994).  Convection with internal sources and turbulence in the earth’s mantle.  Geophys. J. Int, 118:1–19.
  * Tyler, (2008) Tyler, R. H. (2008).  Strong ocean tidal flow and heating on moons of the outer planets.  Nature, 456(7223):770–772.
  * Tyler, (2020) Tyler, R. H. (2020).  Heating of enceladus due to the dissipation of ocean tides.  Icarus, 348:113821.
  * Tyler et al., (2015) Tyler, R. H., Henning, W. G., and Hamilton, C. W. (2015).  Tidal heating in a magma ocean within jupiter’s moon io.  The Astrophysical Journal Supplement Series, 218(2):22.
  * Vance et al., (2007) Vance, S., Harnmeijer, J., Kimura, J., Hussmann, H., DeMartin, B., and Brown, J. M. (2007).  Hydrothermal systems in small ocean planets.  Astrobiology, 7(6):987–1005.
  * Vance et al., (2018) 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 et al., (2020) Vilella, K., Choblet, G., Tsao, W.-E., and Deschamps, F. (2020).  Tidally heated convection and the occurrence of melting in icy satellites: Application to europa.  Journal of Geophysical Research: Planets, 125(3):e2019JE006248.
  * Waite Jr et al., (2009) 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 and Boggs, (2015) Williams, J. G. and Boggs, D. H. (2015).  Tides on the moon: Theory and determination of dissipation.  Journal of Geophysical Research: Planets, 120(4):689–724.

##  Appendix A Sundberg-Cooper viscoelastic model

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!,$
(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)}\bigg{[}\big{(}[J_{U}\eta\omega]^{\alpha}C+1\big{)}\lambda+\frac{J_{U}}{J_{R}}\bigg{]},$
(20)
$\Im[G(\omega)]=\frac{\eta\omega}{H(\omega)}\bigg{[}(J_{U}\eta\omega)^{-1-\alpha}S\lambda+\frac{\eta_{1}}{\eta}+\big{(}\frac{\eta_{1}}{\eta}\big{)}^{2}+(\eta\omega
J_{R})^{-2}\bigg{]},$ (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}}+\bigg{(}\frac{\eta_{1}}{\eta}+1\bigg{)}^{2}+\lambda(J_{U}\eta\omega)^{-2\alpha}(a!)^{2}+2(J_{U}\eta\omega)\bigg{[}S\bigg{(}J_{U}\eta_{1}\omega\big{(}1+\frac{\eta_{1}}{\eta}\big{)}+\frac{J_{U}}{J_{R}(J_{R}\eta\omega)}\bigg{)}+C\big{(}\lambda+\frac{J_{U}}{J_{R}}\big{)}\bigg{]}.$
(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)]}.$ (23)

Parameter (Unit) | Symbol | Value/expression
---|---|---
Constant | S | $\alpha!\sin(\alpha\pi/2)$
Constant | C | $\alpha!\cos(\alpha\pi/2)$
Burgers Coefficient | $\lambda$ | $(J_{U}\eta_{1}\omega)^{2}+(J_{U}/J_{R})^{2}$
Relaxed compliance $\rm(1/Pa.s)$ | $J_{R}$ | $0.2J_{U}$
Voigt-Kelvin Viscosity $\rm(Pa.s)$ | $\eta_{1}$ | $0.02\eta$
Table 4: Parameters associated with the Sundberg-Cooper viscoelastic model.

##  Appendix B Orbital evolution theory

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

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

where $x$ is either $a$ or $e$. The two terms refer to the tides in the planet
and the tides in the satellite, respectively. Because of the spherical shape
of both Pluto and Charon, we can omit the contribution from libration. In what
follows, we provide only the main formulae for the orbital evolution.

The expression for the semi-major axis rate is

$\displaystyle\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)$
(25)
$\displaystyle\bigg{[}\bigg{(}\frac{R}{a}\bigg{)}^{2l+1}\frac{M^{\prime}}{M}F_{lmp}^{2}(i)K_{l}(\beta)+\bigg{(}\frac{R^{\prime}}{a}\bigg{)}^{2l+1}\frac{M}{M^{\prime}}F_{lmp}^{2}(i^{\prime})K_{l}(\beta^{\prime})\bigg{]},$

where $M$ and $M^{\prime}$ are the planet and satellite masses, $R$ and
$R^{{}^{\prime}}$ are their radii, $K_{l}=k_{l}/Q_{l}$ and
$K_{l}^{\,\prime}=k_{l}^{\,\prime}/Q_{l}^{\,\prime}$ are the planet and
satellite quality functions, $G(e)$ and $F(i)$ are the eccentricity and
inclination functions, respectively, $\beta$ 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 are equal to $\frac{1}{2}$ and $3$, respectively.
Similarly to $da/dt$, the general expression for the eccentricity rate is

$\displaystyle\frac{de}{dt}=-\frac{1-e^{2}}{e}\frac{n}{MM^{\prime}}\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)$
(26)
$\displaystyle\bigg{[}\bigg{(}\frac{R}{a}\bigg{)}^{2l+1}\frac{M^{\prime}}{M}F_{lmp}^{2}(i)K_{l}(\beta)+\bigg{(}\frac{R^{\prime}}{a}\bigg{)}^{2l+1}\frac{M}{M^{\prime}}F_{lmp}^{2}(i^{\prime})K^{\prime}_{l}(\beta^{\prime})\bigg{]}+$
$\displaystyle-\frac{\sqrt{1-e^{2}}}{e}\frac{n}{MM^{\prime}}\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)$
$\displaystyle\bigg{[}\bigg{(}\frac{R}{a}\bigg{)}^{2l+1}\frac{M^{\prime}}{M}F_{lmp}^{2}(i)K_{l}(\beta)+\bigg{(}\frac{R^{\prime}}{a}\bigg{)}^{2l+1}\frac{M}{M^{\prime}}F_{lmp}^{2}(i^{\prime})K^{\prime}_{l}(\beta^{\prime})\bigg{]}.$

Similarly, equations for the tidal heating and the tidal torque are given as:

$\displaystyle\bigg{(}\frac{dE}{dt}\bigg{)}_{\text{primary}}=\frac{GM^{\prime
2}}{a}\sum_{l=2}^{\infty}\bigg{(}\frac{R}{a}\bigg{)}^{2l+1}\sum_{m=0}^{l}\frac{(l-m)!}{(l+m)!}(2-\delta_{0m})\sum_{p=0}^{l}\sum_{q=-\infty}^{\infty}G_{lpq}^{2}(e)F_{lmp}^{2}(i)K_{l}(\beta^{\prime})\beta^{\prime},$
(27)

$\displaystyle\bigg{(}\frac{d\dot{\theta}}{dt}\bigg{)}_{\text{primary}}=\frac{GM^{\prime
2}}{Ca}\sum_{l=2}^{\infty}\bigg{(}\frac{R}{a}\bigg{)}^{2l+1}\sum_{m=0}^{l}\frac{(l-m)!}{(l+m)!}(2-\delta_{0m})\sum_{p=0}^{l}\sum_{q=-\infty}^{\infty}G_{lpq}^{2}(e)F_{lmp}^{2}(i)mK_{l}(\beta^{\prime}),$
(28)

where $C$ is the polar moment of inertia of the planet. Similar expressions
are obtained for the satellite are obtained by replacing $M$ with $M^{\prime}$
and $C$ by $C^{\prime}$, where $C^{\prime}$ is the polar moment of inertia of
the secondary.

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

$\displaystyle\bigg{(}\frac{da}{dt}\bigg{)}_{\text{primary}}=n\bigg{(}\frac{R^{5}}{a^{4}}\bigg{)}\bigg{(}\frac{M^{\,\prime}}{M}\bigg{)}\times{\cal{F}}(K_{l},\,\dot{\theta},\,n,\,e)~{}~{}.$
(29)

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

$\bigg{(}\frac{da}{dt}\bigg{)}_{\text{secondary}}=n\bigg{(}\frac{{R^{\,\prime}}^{5}}{a^{4}}\bigg{)}\bigg{(}\frac{M}{M^{\,\prime}}\bigg{)}\times{\cal{F}}(K^{\,\prime}_{l},\,\dot{\theta}^{\,\prime},\,n,\,e)~{}~{}.$
(30)

Here $\,{\cal{F}}\,$ is a function of the eccentricity, the mean motion, the
satellite spin rate ($\dot{\theta}^{\prime}$) and its tidal response

${\cal{F}}\left(K_{l}^{\prime},\dot{\theta}^{\prime},n,e\right)=\sum\limits_{i=0}^{9}e^{2i}\left(\sum\limits_{j=-7}^{-1}K_{l}^{\prime}\left(jn-2\dot{\theta}^{\prime}\right)\varphi_{j}^{2i}+\sum\limits_{j=1}^{11}K_{l}^{\prime}\left(jn-2\dot{\theta}^{\prime}\right)\varphi_{j}^{2i}+\sum\limits_{j=1}^{9}K_{l}^{\prime}(jn)\hat{\varphi}_{j}^{2i}\right).$
(31)

The coefficients $\varphi_{j}^{2i}$ and $\hat{\varphi}_{j}^{2i}$ of the series
are tabulated in Supplementary Tables 3 and 4. 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 synchronised
moon, we have $\,\dot{\theta}^{\,\prime}=n\,$, wherefore the semi-diurnal term
in equation 30 vanishes: $\,K_{l}(2n-2\dot{\theta}^{\,\prime})\,=\,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}^{\prime}\,$.

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

$\displaystyle\bigg{(}\frac{de}{dt}\bigg{)}_{\text{primary}}=-n\frac{M^{\prime}}{M}\bigg{(}\frac{R}{a}\bigg{)}^{5}\times{\cal{L}}(K_{l},\dot{\theta},n,e)~{}~{},$
(32)

while the input from the tides in the satellite reads as

$\displaystyle\bigg{(}\frac{de}{dt}\bigg{)}_{\text{secondary}}=-n\frac{M}{M^{\prime}}\bigg{(}\frac{R^{\prime}}{a}\bigg{)}^{5}\times{\cal{L}}(K^{\,\prime}_{l},\dot{\theta}^{\,\prime},n,e)~{}~{},$
(33)

where the function $\cal{L}$ is defined as

${\cal
L}(K_{l},\dot{\theta},n,e)=\sum\limits_{i=1}^{9}e^{2i-1}\left(\sum\limits_{j=-7}^{-1}K_{l}\left(jn-2\dot{\theta}\right)\lambda_{j}^{2i-1}+\sum\limits_{j=1}^{11}K_{l}\left(jn-2\dot{\theta}\right)\lambda_{j}^{2i-1}+\sum\limits_{j=1}^{9}K_{l}(jn)\hat{\lambda}_{j}^{2i-1}\right)~{}~{}.$
(34)

The coefficients $\lambda_{j}^{2i-1}$ and $\hat{\lambda}_{j}^{2i-1}$ are
tabulated in Supplementary Table 6 and 7. Note that, similarly to $da/dt$, the
expression is general, in that neither of the bodies is a priori assumed
synchronous. For the synchronised case, the term associated with the semi-
diurnal tide in equation (32) vanishes.

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

$\displaystyle\bigg{(}\frac{dE}{dt}\bigg{)}_{\text{primary}}=\frac{GM^{\prime
2}}{a}\bigg{(}\frac{R}{a}\bigg{)}^{5}\times{\cal{H}}(K_{l},\dot{\theta},n,e)~{}~{},$
(35)

The expression for ${\cal{H}}$ is given by

${\cal{H}}(K_{l},\dot{\theta},n,e)=\sum\limits_{i=0}^{9}e^{2i}\left[\sum\limits_{j=-7}^{-1}\xi_{j}^{2i}|jn-2\dot{\theta}|K_{l}\left(|jn-2\dot{\theta}|\right)+\sum\limits_{j=1}^{11}\xi_{j}^{2i}|jn-2\dot{\theta}|K_{l}\left(|jn-2\dot{\theta}|\right)+\sum\limits_{j=1}^{9}jn\hat{\xi}_{j}^{2i}K_{l}(jn)\right]$
(36)

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

$\displaystyle\bigg{(}\frac{d\dot{\theta}}{dt}\bigg{)}_{\text{primary}}=\frac{GM^{\prime}}{Ca}\bigg{(}\frac{R}{a}\bigg{)}^{5}\times{\cal{G}}(K_{l},\dot{\theta},n,e)~{}~{},$
(37)

The expression ${\cal{G}}$ is

${\cal{G}}(K_{l},\dot{\theta},n,e)=\sum\limits_{i=0}^{9}e^{2i}\left[\sum\limits_{j=-7}^{-1}\chi_{j}^{2i}K_{l}\left(jn-2\dot{\theta}\right)+\sum\limits_{j=1}^{11}\chi_{j}^{2i}K_{l}\left(jn-2\dot{\theta}\right)\right]$
(38)

## Supplementary material: Tidal coefficients

This section summarizes all the tidal coefficients that are required in the
development of the orbital theory as formulated in the Methods section of the
main text, i.e., Eq. (15)–(18).

Table 5: $\varphi_{j}^{2i}$ $\varphi_{j}^{2i}$
---
${}_{j}\backslash^{2i}$ | $0$ | $2$ | $4$ | $6$
$-1$ | $0$ | $0$ | $0$ | $\frac{1}{1536}$
$1$ | $0$ | $-\frac{3}{8}$ | $\frac{3}{32}$ | $-\frac{13}{512}$
$2$ | $-3$ | $15$ | $-\frac{189}{8}$ | $\frac{155}{12}$
$3$ | $0$ | $-\frac{441}{8}$ | $\frac{7749}{32}$ | $-\frac{197775}{512}$
$4$ | $0$ | $0$ | $-\frac{867}{2}$ | $1955$
$5$ | $0$ | $0$ | $0$ | $-\frac{3570125}{1536}$

$\varphi_{j}^{2i}$
---
${}_{j}\backslash{2i}$ | $8$ | $10$ | $12$ | $14$ | $16$ | $18$
$-7$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{232630513987207}{469654673817600}$
$-6$ | $0$ | $0$ | $0$ | $0$ | $\frac{4782969}{20070400}$ | $\frac{1594323}{20070400}$
$-5$ | $0$ | $0$ | $0$ | $\frac{1220703125}{11098128384}$ | $\frac{13427734375}{177570054144}$ | $\frac{1619873046875}{14611478740992}$
$-4$ | $0$ | $0$ | $\frac{32}{675}$ | $\frac{32}{675}$ | $\frac{92}{1575}$ | $\frac{7808}{127575}$
$-3$ | $0$ | $\frac{59049}{3276800}$ | $\frac{59049}{2621440}$ | $\frac{38558997}{1468006400}$ | $\frac{18482337}{671088640}$ | $\frac{72655602021}{2630667468800}$
$-2$ | $\frac{1}{192}$ | $\frac{7}{960}$ | $\frac{949}{115200}$ | $\frac{10193}{1209600}$ | $\frac{127453}{15482880}$ | $\frac{781813}{99532800}$
$-1$ | $\frac{11}{12288}$ | $\frac{619}{655360}$ | $\frac{64537}{70778880}$ | $\frac{33683687}{39636172800}$ | $\frac{165074183}{211392921600}$ | $\frac{195618433037}{273965226393600}$
$1$ | $-\frac{113}{12288}$ | $-\frac{2639}{327680}$ | $-\frac{129709}{19660800}$ | $-\frac{73185503}{13212057600}$ | $-\frac{21004537637}{4439251353600}$ | $-\frac{5242619992367}{1278504389836800}$
$2$ | $-\frac{2881}{768}$ | $\frac{3481}{6400}$ | $-\frac{12791}{92160}$ | $-\frac{288163}{5644800}$ | $-\frac{1190849}{24084480}$ | $-\frac{535500317}{12541132800}$
$3$ | $\frac{1193715}{4096}$ | $-\frac{41889501}{327680}$ | $\frac{237091671}{6553600}$ | $-\frac{1586248587}{209715200}$ | $\frac{863843049}{939524096}$ | $-\frac{220413818169}{751619276800}$
$4$ | $-\frac{83551}{24}$ | $\frac{131319}{40}$ | $-\frac{2459537}{1280}$ | $\frac{232033433}{302400}$ | $-\frac{77502109}{345600}$ | $\frac{1505096129}{30481920}$
$5$ | $\frac{137418125}{12288}$ | $-\frac{2936126875}{131072}$ | $\frac{356500341875}{14155776}$ | $-\frac{29029068521875}{1585446912}$ | $\frac{79231092391625}{8455716864}$ | $-\frac{5599532275549375}{1565515579392}$
$6$ | $-\frac{2556801}{256}$ | $\frac{66328119}{1280}$ | $-\frac{5938832817}{51200}$ | $\frac{13527027819}{89600}$ | $-\frac{1497798429597}{11468800}$ | $\frac{931614290433}{11468800}$
$7$ | $0$ | $-\frac{364996466863}{9830400}$ | $\frac{981779068235}{4718592}$ | $-\frac{2921128205246293}{5662310400}$ | $\frac{2760195487312025}{3623878656}$ | $-\frac{9940709866552998323}{13045963161600}$
$8$ | $0$ | $0$ | $-\frac{5383010161}{43200}$ | $\frac{56914314263}{75600}$ | $-\frac{578658802849}{282240}$ | $\frac{77318558814617}{22861440}$
$9$ | $0$ | $0$ | $0$ | $-\frac{3982050900860283}{10276044800}$ | $\frac{82748209967119359}{32883343360}$ | $-\frac{11230592195266666389}{1503238553600}$
$10$ | $0$ | $0$ | $0$ | $0$ | $-\frac{2010318936125}{1769472}$ | $\frac{880939918004375}{111476736}$
$11$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{73137037043155488990059}{23013079017062400}$

Table 6: $\hat{\varphi}_{j}^{2i}$ $\hat{\varphi}_{j}^{2i}$
---
${}_{j}\backslash^{2i}$ | $0$ | $2$ | $4$ | $6$
$1$ | $0$ | $-\frac{9}{4}$ | $-\frac{81}{16}$ | $-\frac{2295}{256}$
$2$ | $0$ | $0$ | $-\frac{81}{8}$ | $-\frac{63}{4}$
$3$ | $0$ | $0$ | $0$ | $-\frac{8427}{256}$
$\hat{\varphi}_{j}^{2i}$
---
${}_{j}\backslash^{2i}$ | $8$ | $10$ | $12$ | $14$ | $16$ | $18$
$1$ | $-\frac{28403}{2048}$ | $-\frac{3240741}{163840}$ | $-\frac{87480539}{3276800}$ | $-\frac{76230921863}{2202009600}$ | $-\frac{143189681823}{3288334336}$ | $-\frac{3798168072196597}{71028021657600}$
$2$ | $-\frac{1661}{64}$ | $-\frac{12027}{320}$ | $-\frac{525521}{10240}$ | $-\frac{13519469}{201600}$ | $-\frac{729537581}{8601600}$ | $-\frac{56670971611}{541900800}$
$3$ | $-\frac{62487}{2048}$ | $-\frac{18059643}{327680}$ | $-\frac{96458913}{1310720}$ | $-\frac{10219795947}{104857600}$ | $-\frac{1455222699297}{11744051200}$ | $-\frac{28832512267869}{187904819200}$
$4$ | $-\frac{5929}{64}$ | $-\frac{9933}{320}$ | $-\frac{708871}{6400}$ | $-\frac{1736021}{14400}$ | $-\frac{27865091}{172032}$ | $-\frac{2578148063}{12902400}$
$5$ | $0$ | $-\frac{15717645}{65536}$ | $\frac{14736585}{262144}$ | $-\frac{67620980465}{264241152}$ | $-\frac{212874343975}{1409286144}$ | $-\frac{13593131130225}{52613349376}$
$6$ | $0$ | $0$ | $-\frac{30089667}{51200}$ | $\frac{20778687}{44800}$ | $-\frac{592612065}{802816}$ | $-\frac{124700799}{20070400}$
$7$ | $0$ | $0$ | $0$ | $-\frac{1306918425967}{943718400}$ | $\frac{9297259269991}{5033164800}$ | $-\frac{6999935835286093}{2899102924800}$
$8$ | $0$ | $0$ | $0$ | $0$ | $-\frac{31801945561}{10035200}$ | $\frac{1606626956771}{270950400}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{37292291277118641}{5261334937600}$
Table 7: $\gamma_{0}^{j,s}$ $\gamma_{0}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$ | $10$ | $12$ | $14$ | $16$ | $18$
$1$ | $-3$ | $3$ | $-\frac{5}{4}$ | $\frac{7}{24}$ | $-\frac{7}{160}$ | $\frac{11}{2400}$ | $-\frac{143}{403200}$ | $\frac{143}{6773760}$ | $-\frac{2431}{2438553600}$
$2$ | $0$ | $-\frac{3}{2}$ | $1$ | $-\frac{7}{24}$ | $\frac{1}{20}$ | $-\frac{11}{1920}$ | $\frac{143}{302400}$ | $-\frac{143}{4838400}$ | $\frac{221}{152409600}$
$3$ | $0$ | $0$ | $-\frac{1}{4}$ | $\frac{1}{8}$ | $-\frac{9}{320}$ | $\frac{11}{2880}$ | $-\frac{143}{403200}$ | $\frac{13}{537600}$ | $-\frac{221}{174182400}$
$4$ | $0$ | $0$ | $0$ | $-\frac{1}{48}$ | $\frac{1}{120}$ | $-\frac{11}{7200}$ | $\frac{13}{75600}$ | $-\frac{13}{967680}$ | $\frac{17}{21772800}$
$5$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1}{960}$ | $\frac{1}{2880}$ | $-\frac{13}{241920}$ | $\frac{1}{193536}$ | $-\frac{17}{48771072}$
$6$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1}{28800}$ | $\frac{1}{100800}$ | $-\frac{1}{752640}$ | $\frac{17}{152409600}$
$7$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1}{1209600}$ | $\frac{1}{4838400}$ | $-\frac{17}{696729600}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1}{67737600}$ | $\frac{1}{304819200}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1}{4877107200}$

$\gamma_{1}^{j,s}$
---
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$ | $9$ | $11$ | $13$ | $15$ | $17$
$1$ | $30$ | $-42$ | $\frac{45}{2}$ | $-\frac{77}{12}$ | $\frac{91}{80}$ | $-\frac{11}{80}$ | $\frac{2431}{201600}$ | $-\frac{2717}{3386880}$ | $\frac{2431}{58060800}$
$2$ | $0$ | $21$ | $-18$ | $\frac{77}{12}$ | $-\frac{13}{10}$ | $\frac{11}{64}$ | $-\frac{2431}{151200}$ | $\frac{2717}{2419200}$ | $-\frac{221}{3628800}$
$3$ | $0$ | $0$ | $\frac{9}{2}$ | $-\frac{11}{4}$ | $\frac{117}{160}$ | $-\frac{11}{96}$ | $\frac{2431}{201600}$ | $-\frac{247}{268800}$ | $\frac{221}{4147200}$
$4$ | $0$ | $0$ | $0$ | $\frac{11}{24}$ | $-\frac{13}{60}$ | $\frac{11}{240}$ | $-\frac{221}{37800}$ | $\frac{247}{483840}$ | $-\frac{17}{518400}$
$5$ | $0$ | $0$ | $0$ | $0$ | $\frac{13}{480}$ | $-\frac{1}{96}$ | $\frac{221}{120960}$ | $-\frac{19}{96768}$ | $\frac{17}{1161216}$
$6$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{1}{960}$ | $-\frac{17}{50400}$ | $\frac{19}{376320}$ | $-\frac{17}{3628800}$
$7$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{17}{604800}$ | $-\frac{19}{2419200}$ | $\frac{17}{16588800}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{19}{33868800}$ | $-\frac{1}{7257600}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{1}{116121600}$

$\gamma_{2}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$ | $10$ | $12$ | $14$ | $16$
$1$ | $\frac{825}{4}$ | $-\frac{2753}{16}$ | $\frac{6313}{96}$ | $-\frac{9303}{640}$ | $\frac{20153}{9600}$ | $-\frac{68959}{322560}$ | $\frac{2192333}{135475200}$ | $-\frac{9231079}{9754214400}$
$2$ | $-\frac{249}{2}$ | $\frac{305}{2}$ | $-\frac{421}{6}$ | $\frac{139}{8}$ | $-\frac{5207}{1920}$ | $\frac{176803}{604800}$ | $-\frac{55913}{2419200}$ | $\frac{426439}{304819200}$
$3$ | $0$ | $-\frac{709}{16}$ | $\frac{1063}{32}$ | $-\frac{2685}{256}$ | $\frac{21953}{11520}$ | $-\frac{52613}{230400}$ | $\frac{14001}{716800}$ | $-\frac{125099}{99532800}$
$4$ | $0$ | $0$ | $-\frac{301}{48}$ | $\frac{817}{240}$ | $-\frac{2941}{3600}$ | $\frac{101}{864}$ | $-\frac{54847}{4838400}$ | $\frac{4987}{6220800}$
$5$ | $0$ | $0$ | $0$ | $-\frac{1817}{3840}$ | $\frac{2323}{11520}$ | $-\frac{37753}{967680}$ | $\frac{17803}{3870720}$ | $-\frac{72977}{195084288}$
$6$ | $0$ | $0$ | $0$ | $0$ | $-\frac{637}{28800}$ | $\frac{313}{40320}$ | $-\frac{789}{627200}$ | $\frac{38419}{304819200}$
$7$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{3397}{4838400}$ | $\frac{811}{3870720}$ | $-\frac{81401}{2786918400}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1091}{67737600}$ | $\frac{2549}{609638400}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{5449}{19508428800}$
$\gamma_{3}^{j,s}$
---
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$ | $9$ | $11$ | $13$ | $15$
$1$ | $-\frac{1647}{4}$ | $\frac{5561}{8}$ | $-\frac{38129}{96}$ | $\frac{11165}{96}$ | $-\frac{40039}{1920}$ | $\frac{2042161}{806400}$ | $-\frac{4303871}{19353600}$ | $\frac{6010433}{406425600}$
$2$ | $357$ | $-\frac{2993}{4}$ | $\frac{23023}{48}$ | $-\frac{36377}{240}$ | $\frac{11029}{384}$ | $-\frac{4400891}{1209600}$ | $\frac{639353}{1935360}$ | $-\frac{655343}{29030400}$
$3$ | $0$ | $\frac{2137}{8}$ | $-\frac{8463}{32}$ | $\frac{32879}{320}$ | $-\frac{50911}{2304}$ | $\frac{2464253}{806400}$ | $-\frac{1911013}{6451200}$ | $\frac{11453}{537600}$
$4$ | $0$ | $0$ | $\frac{5551}{96}$ | $-\frac{3619}{96}$ | $\frac{7553}{720}$ | $-\frac{515197}{302400}$ | $\frac{595907}{3225600}$ | $-\frac{419393}{29030400}$
$5$ | $0$ | $0$ | $0$ | $\frac{1881}{320}$ | $-\frac{6605}{2304}$ | $\frac{300767}{483840}$ | $-\frac{104549}{1290240}$ | $\frac{2095}{290304}$
$6$ | $0$ | $0$ | $0$ | $0$ | $\frac{1991}{5760}$ | $-\frac{54389}{403200}$ | $\frac{77323}{3225600}$ | $-\frac{19639}{7526400}$
$7$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{31991}{2419200}$ | $-\frac{83323}{19353600}$ | $\frac{75289}{116121600}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{6871}{19353600}$ | $-\frac{40321}{406425600}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{13}{1843200}$
$\gamma_{4}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$ | $10$ | $12$ | $14$
$1$ | $-\frac{47555}{32}$ | $\frac{190709}{128}$ | $-\frac{400233}{640}$ | $\frac{1121033}{7680}$ | $-\frac{23458073}{1075200}$ | $\frac{35199901}{15482880}$ | $-\frac{852863297}{4877107200}$
$2$ | $\frac{8687}{4}$ | $-\frac{34773}{16}$ | $\frac{4644}{5}$ | $-\frac{425111}{1920}$ | $\frac{1512829}{44800}$ | $-\frac{1387199}{387072}$ | $\frac{85408297}{304819200}$
$3$ | $-\frac{32127}{32}$ | $\frac{185441}{128}$ | $-\frac{934671}{1280}$ | $\frac{8842181}{46080}$ | $-\frac{100865411}{3225600}$ | $\frac{9994067}{2867200}$ | $-\frac{98657507}{348364800}$
$4$ | $0$ | $-\frac{23833}{64}$ | $\frac{12277}{40}$ | $-\frac{196591}{1920}$ | $\frac{1944769}{100800}$ | $-\frac{9191771}{3870720}$ | $\frac{9035869}{43545600}$
$5$ | $0$ | $0$ | $-\frac{68941}{1280}$ | $\frac{31835}{1024}$ | $-\frac{998839}{129024}$ | $\frac{17575681}{15482880}$ | $-\frac{10930351}{97542144}$
$6$ | $0$ | $0$ | $0$ | $-\frac{23699}{5760}$ | $\frac{737969}{403200}$ | $-\frac{5213}{14336}$ | $\frac{13298189}{304819200}$
$7$ | $0$ | $0$ | $0$ | $0$ | $-\frac{208333}{1075200}$ | $\frac{5434133}{77414400}$ | $-\frac{16170799}{1393459200}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{239213}{38707200}$ | $\frac{1152077}{609638400}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{1387327}{9754214400}$
$\gamma_{5}^{j,s}$
---
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$ | $9$ | $11$ | $13$
$1$ | $\frac{22769}{16}$ | $-\frac{1296283}{384}$ | $\frac{4380193}{1920}$ | $-\frac{42520013}{57600}$ | $\frac{113941067}{806400}$ | $-\frac{463361041}{25804800}$ | $\frac{1135527679}{696729600}$
$2$ | $-\frac{13539}{4}$ | $\frac{2484625}{384}$ | $-\frac{1952053}{480}$ | $\frac{58912711}{46080}$ | $-\frac{292763719}{1209600}$ | $\frac{2374678163}{77414400}$ | $-\frac{53957761}{19353600}$
$3$ | $\frac{71825}{32}$ | $-\frac{697647}{128}$ | $\frac{9779869}{2560}$ | $-\frac{14726891}{11520}$ | $\frac{813614173}{3225600}$ | $-\frac{850321511}{25804800}$ | $\frac{854790599}{278691840}$
$4$ | $0$ | $\frac{81433}{48}$ | $-\frac{3590417}{1920}$ | $\frac{89371417}{115200}$ | $-\frac{210499391}{1209600}$ | $\frac{958855703}{38707200}$ | $-\frac{854451569}{348364800}$
$5$ | $0$ | $0$ | $\frac{191637}{512}$ | $-\frac{303739}{1152}$ | $\frac{149339629}{1935360}$ | $-\frac{200654617}{15482880}$ | $\frac{400813631}{278691840}$
$6$ | $0$ | $0$ | $0$ | $\frac{8889011}{230400}$ | $-\frac{8043643}{403200}$ | $\frac{116210893}{25804800}$ | $-\frac{105142393}{174182400}$
$7$ | $0$ | $0$ | $0$ | $0$ | $\frac{11083823}{4838400}$ | $-\frac{10353229}{11059200}$ | $\frac{159751717}{928972800}$
$8$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{2277181}{25804800}$ | $-\frac{20783387}{696729600}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{1327271}{557383680}$

$\gamma_{6}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$ | $10$ | $12$
$1$ | $\frac{6360697}{1536}$ | $-\frac{83914747}{15360}$ | $\frac{205972691}{76800}$ | $-\frac{238333495}{344064}$ | $\frac{34401989173}{309657600}$ | $-\frac{67943553227}{5573836800}$
$2$ | $-\frac{4527427}{384}$ | $\frac{11995117}{960}$ | $-\frac{85032277}{15360}$ | $\frac{218343289}{161280}$ | $-\frac{1086174427}{5160960}$ | $\frac{7897450723}{348364800}$
$3$ | $\frac{6894791}{512}$ | $-\frac{149105731}{10240}$ | $\frac{67081587}{10240}$ | $-\frac{2785868387}{1720320}$ | $\frac{2911274687}{11468800}$ | $-\frac{25580115979}{928972800}$
$4$ | $-\frac{1024813}{192}$ | $\frac{818159}{96}$ | $-\frac{22057967}{4800}$ | $\frac{1273775}{1008}$ | $-\frac{4112465311}{19353600}$ | $\frac{421956389}{17418240}$
$5$ | $0$ | $-\frac{20343329}{10240}$ | $\frac{54424859}{30720}$ | $-\frac{644732521}{1032192}$ | $\frac{7580219749}{61931520}$ | $-\frac{8635714091}{557383680}$
$6$ | $0$ | $0$ | $-\frac{7421471}{25600}$ | $\frac{1913875}{10752}$ | $-\frac{133004803}{2867200}$ | $\frac{816994861}{116121600}$
$7$ | $0$ | $0$ | $0$ | $-\frac{115244509}{5160960}$ | $\frac{1077524291}{103219200}$ | $-\frac{960692647}{445906944}$
$8$ | $0$ | $0$ | $0$ | $0$ | $-\frac{5122889}{4838400}$ | $\frac{1242901}{3110400}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $-\frac{18025373}{530841600}$
$\gamma_{7}^{j,s}$
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$ | $9$ | $11$
$1$ | $-\frac{1663559}{768}$ | $\frac{121709581}{15360}$ | $-\frac{312427823}{46080}$ | $\frac{19673089399}{7741440}$ | $-\frac{291425548999}{541900800}$ | $\frac{2867948486117}{39016857600}$
$2$ | $\frac{4244105}{384}$ | $-\frac{195851269}{7680}$ | $\frac{811588697}{46080}$ | $-\frac{12636131167}{2150400}$ | $\frac{1402888481}{1209600}$ | $-\frac{986035976887}{6502809600}$
$3$ | $-\frac{9859807}{512}$ | $\frac{24914363}{640}$ | $-\frac{4694903647}{184320}$ | $\frac{11825097051}{1433600}$ | $-\frac{23555895221}{14745600}$ | $\frac{1149280143791}{5573836800}$
$4$ | $\frac{4053463}{384}$ | $-\frac{3400213}{120}$ | $\frac{1948423099}{92160}$ | $-\frac{317943461}{43008}$ | $\frac{11099163623}{7372800}$ | $-\frac{280439920391}{1393459200}$
$5$ | $0$ | $\frac{121754083}{15360}$ | $-\frac{1736132939}{184320}$ | $\frac{31935948797}{7741440}$ | $-\frac{19885290673}{20643840}$ | $\frac{1102592661521}{7803371520}$
$6$ | $0$ | $0$ | $\frac{8062469}{4608}$ | $-\frac{1689833129}{1290240}$ | $\frac{4530651799}{11289600}$ | $-\frac{1362209431217}{19508428800}$
$7$ | $0$ | $0$ | $0$ | $\frac{7010016683}{38707200}$ | $-\frac{15241243613}{154828800}$ | $\frac{14310329239}{619315200}$
$8$ | $0$ | $0$ | $0$ | $0$ | $\frac{23423982107}{2167603200}$ | $-\frac{89822032001}{19508428800}$
$9$ | $0$ | $0$ | $0$ | $0$ | $0$ | $\frac{2039795207}{4877107200}$

$\gamma_{8}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$ | $10$
$1$ | $-\frac{36667771}{6144}$ | $\frac{137576377}{12288}$ | $-\frac{58268412071}{8601600}$ | $\frac{3483804392591}{1734082560}$ | $-\frac{3962944769797}{11147673600}$
$2$ | $\frac{7932001}{256}$ | $-\frac{3629434021}{92160}$ | $\frac{62558288021}{3225600}$ | $-\frac{393958246799}{77414400}$ | $\frac{1159836156341}{1393459200}$
$3$ | $-\frac{1378786789}{20480}$ | $\frac{5365102241}{73728}$ | $-\frac{847780261361}{25804800}$ | $\frac{3370548827407}{412876800}$ | $-\frac{5974574495}{4644864}$
$4$ | $\frac{124022959}{1920}$ | $-\frac{425217367}{5760}$ | $\frac{13920972659}{403200}$ | $-\frac{34003821749}{3870720}$ | $\frac{122575184351}{87091200}$
$5$ | $-\frac{1391241377}{61440}$ | $\frac{4788563821}{122880}$ | $-\frac{114133290857}{5160960}$ | $\frac{1565177035693}{247726080}$ | $-\frac{609862623793}{557383680}$
$6$ | $0$ | $-\frac{154509701}{18432}$ | $\frac{3653368319}{460800}$ | $-\frac{105576999331}{36126720}$ | $\frac{275386936571}{464486400}$
$7$ | $0$ | $0$ | $-\frac{31536438929}{25804800}$ | $\frac{975977916433}{1238630400}$ | $-\frac{4758029916799}{22295347200}$
$8$ | $0$ | $0$ | $0$ | $-\frac{25549718947}{270950400}$ | $\frac{4573097093}{99532800}$
$9$ | $0$ | $0$ | $0$ | $0$ | $-\frac{99964117517}{22295347200}$
$\gamma_{9}^{j,s}$
---
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$ | $9$
$1$ | $\frac{72574433}{40960}$ | $-\frac{33250133}{3072}$ | $\frac{321566834083}{25804800}$ | $-\frac{48892934193757}{8670412800}$ | $\frac{424304933869777}{312134860800}$
$2$ | $-\frac{139238303}{7680}$ | $\frac{2812070945}{49152}$ | $-\frac{1815292535809}{38707200}$ | $\frac{21601089607747}{1238630400}$ | $-\frac{144839901437071}{39016857600}$
$3$ | $\frac{326578731}{5120}$ | $-\frac{3318429817}{23040}$ | $\frac{1712182537183}{17203200}$ | $-\frac{13787028701293}{412876800}$ | $\frac{295876625683021}{44590694400}$
$4$ | $-\frac{45278785}{512}$ | $\frac{8612926109}{46080}$ | $-\frac{4880466678209}{38707200}$ | $\frac{17258868186769}{412876800}$ | $-\frac{22984018946099}{2786918400}$
$5$ | $\frac{851962171}{20480}$ | $-\frac{44593891643}{368640}$ | $\frac{293291537465}{3096576}$ | $-\frac{8497199910629}{247726080}$ | $\frac{89611617850873}{12485394432}$
$6$ | $0$ | $\frac{713014643}{23040}$ | $-\frac{41944969211}{1075200}$ | $\frac{51535195929577}{2890137600}$ | $-\frac{167924593843033}{39016857600}$
$7$ | $0$ | $0$ | $\frac{1052720951633}{154828800}$ | $-\frac{6613213406581}{1238630400}$ | $\frac{151464161578421}{89181388800}$
$8$ | $0$ | $0$ | $0$ | $\frac{76155429263}{108380160}$ | $-\frac{31037807709701}{78033715200}$
$9$ | $0$ | $0$ | $0$ | $0$ | $\frac{26184666395077}{624269721600}$
$\gamma_{10}^{j,s}$
---
${}_{j}\backslash^{s}$ | $2$ | $4$ | $6$ | $8$
$1$ | $\frac{4188185133}{819200}$ | $-\frac{668034376553}{45875200}$ | $\frac{298635437603}{26214400}$ | $-\frac{5008680379354717}{1248539443200}$
$2$ | $-\frac{69622578569}{1474560}$ | $\frac{288361412803}{3686400}$ | $-\frac{3733878109439}{82575360}$ | $\frac{14728473880769}{1114767360}$
$3$ | $\frac{105043267843}{589824}$ | $-\frac{90447675495017}{412876800}$ | $\frac{176246672641693}{1651507200}$ | $-\frac{2489696844203003}{89181388800}$
$4$ | $-\frac{72784945777}{230400}$ | $\frac{560166907039}{1612800}$ | $-\frac{49295772110719}{309657600}$ | $\frac{13961643516553}{348364800}$
$5$ | $\frac{776330530877}{2949120}$ | $-\frac{1732963120055}{5505024}$ | $\frac{30099792470603}{198180864}$ | $-\frac{705923122737949}{17836277760}$
$6$ | $-\frac{613186054801}{7372800}$ | $\frac{436748965809}{2867200}$ | $-\frac{37201323168383}{412876800}$ | $\frac{1036259700960391}{39016857600}$
$7$ | $0$ | $-\frac{4206593801911}{137625600}$ | $\frac{7158257783303}{235929600}$ | $-\frac{4137506209801739}{356725555200}$
$8$ | $0$ | $0$ | $-\frac{229049541811}{51609600}$ | $\frac{58107407588357}{19508428800}$
$9$ | $0$ | $0$ | $0$ | $-\frac{645135681793}{1887436800}$
$\gamma_{11}^{j,s}$
---
${}_{j}\backslash^{s}$ | $1$ | $3$ | $5$ | $7$
$1$ | $-\frac{4390861889}{4915200}$ | $\frac{1961183621879}{206438400}$ | $-\frac{38437694104309}{2477260800}$ | $\frac{5499256218016573}{624269721600}$
$2$ | $\frac{12960937279}{737280}$ | $-\frac{25335942630283}{309657600}$ | $\frac{1196416018969}{14155776}$ | $-\frac{1266083597711297}{34681651200}$