# Thermal Processing of Jupiter-family Comets during Their Chaotic Orbital Evolution

Anastasios Gkotsinas, Aurélie Guilbert-Lepoutre, Sean N. Raymond, David Nesvorny

► To cite this version:

Anastasios Gkotsinas, Aurélie Guilbert-Lepoutre, Sean N. Raymond, David Nesvorny. Thermal Pro- cessing of Jupiter-family Comets during Their Chaotic Orbital Evolution. The Astrophysical Journal, 2022, 928, 10.3847/1538-4357/ac54ac. insu-03678902

HAL Id: insu-03678902

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

Submitted on 25 May 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# Thermal Processing of Jupiter-family Comets during Their Chaotic Orbital Evolution

Anastasios Gkotsinas¹ ID, Aurélie Guilbert-Lepoutre¹ ID, Sean N. Raymond² ID, and David Nesvorny³ ID

¹Laboratoire de Géologie de Lyon: Terre, Planètes, Environnement, CNRS, UCBL, ENSL, F-69622, Villeurbanne, France; anastasios.gkotsinas@univ-lyon1.fr

²Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, F-33615 Pessac, France

³Department of Space Studies, Southwest Research Institute, 1050 Walnut Street, Suite 300, Boulder, CO 80302, USA

Received 2021 December 20; revised 2022 February 8; accepted 2022 February 10; published 2022 March 25

## Abstract

Evidence for cometary activity beyond Jupiter's and Saturn's orbits—such as that observed for Centaurs and long-period comets—suggests that the thermal processing of comet nuclei starts long before they enter the inner solar system, where they are typically observed and monitored. Such observations raise questions as to the depth of unprocessed material and whether the activity of Jupiter-family comets (JFCs) can be representative of any primitive material. Here we model the coupled thermal and dynamical evolution of JFCs, from the moment they leave their outer solar system reservoirs until their ejection into interstellar space. We apply a thermal evolution model to a sample of simulated JFCs obtained from dynamical simulations that successfully reproduce the orbital distribution of observed JFCs. We show that due to the stochastic nature of comet trajectories toward the inner solar system, all simulated JFCs undergo multiple heating episodes resulting in significant modifications of their initial volatile contents. A statistical analysis constrains the extent of such processing. We suggest that primordial condensed hypervolatile ices should be entirely lost from the layers that contribute to cometary activity observed today. Our results demonstrate that understanding the orbital (and thus, heating) history of JFCs is essential when putting observations in a broader context.

Unified Astronomy Thesaurus concepts: Comets (280); Short period comets (1452); Comet nuclei (2160); Comet dynamics (2213); Comet volatiles (2162); Computational methods (1965)

## 1. Introduction

Jupiter-family comets (JFCs) represent a population of icy objects whose orbits are primarily determined by the gravitational influence of Jupiter. JFCs are characterized by short orbital periods (≤ 20 yr) and low inclinations (≤30°) (Di Sisto et al. 2009; Nesvorný et al. 2017). They are thought to originate from the Kuiper Belt and scattered disk (Brasser & Morbidelli 2013) and evolve through the giant-planet region on unstable orbits, before reaching the inner solar system after a close encounter with Jupiter (Levison & Duncan 1997; Di Sisto et al. 2009; Nesvorný et al. 2017; Fernández et al. 2018; Steckloff et al. 2020).

The thermally induced processing of JFCs' physical and chemical properties can be divided into four distinct phases (e.g., Meech & Svoren 2004; Prialnik et al. 2004):

1. An initial stage that encompasses the formation of cometary nuclei themselves and their processing prior to their displacement in the outer solar system reservoirs;

2. A "reservoir" phase, lasting several billion years, where nuclei are thought to remain relatively unaltered;

3. An intermediate phase of orbital perturbations bringing comet nuclei from their reservoirs to the inner solar system;

4. A short-lived phase of intense processing known as the "active phase," where cometary activity is mostly driven by water-ice sublimation at distances within Jupiter's orbit (~3 au).

Although the active phase is considered to drive the most intense and rapid evolution, multiple lines of evidence suggest that accounting for that stage of processing alone does not provide the full picture of the comets' thermal processing. Considerable modeling efforts have examined the possibility of thermal processing and activity even from the earlier stages of evolution. Such studies propose that alterations start as early as the formation of comet nuclei, and before their displacement in outer solar system reservoirs (e.g. Prialnik & Merk 2008; Raymond et al. 2020; Davidsson 2021), but also during the prolonged storage phase, where surface temperatures between 30 and 50 K can be attained (e.g., Capria et al. 2000; De Sanctis et al. 2000, 2001; Choi et al. 2002). The existence of active Centaurs (e.g., Jewitt 2009; Lin et al. 2014; Mazzotta Epifani et al. 2017, 2018; Steckloff et al. 2020; de la Fuente Marcos et al. 2021), alongside the recent observations of the long-distant activity of comets (e.g., Meech et al. 2017; Jewitt et al. 2017; Hui et al. 2018, 2019; Yang et al. 2021; Farnham et al. 2021), provide direct evidence that cometary activity starts beyond the orbit of Neptune.

In this paper, we investigate the thermal processing of JFCs during their complete orbital evolution. We start from the time that cometary nuclei leave their reservoirs and follow their dynamical evolution until they are ejected from the solar system. We apply a thermal evolution model to a sample of model JFCs from the dynamical simulations of Nesvorný et al. (2017). In this sample, model JFCs evolve from Neptune- to Jupiter-crossing orbits with perihelion distances $q < 2.5$ au, reproducing the orbital distribution of currently observed JFCs. We perform thermal evolution calculations coupled to the dynamical evolution of 276 model JFCs in order to obtain the internal temperature distribution and its evolution. This enables us to assess the degree of thermal processing among JFCs in a statistically significant manner.

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In Section 2, we describe our thermal evolution model and our model JFC sample, along with considerations necessary for their coupling. In Section 3, we present the results of our simulations for both individual simulated JFCs and the whole sample. In Section 4, we discuss the limitations of our approach and their consequences on the results, while in Section 5 we present our conclusions.

## 2. Methods

Coupling the thermal and orbital evolution of comet nuclei is not a trivial procedure, as the underlying processes act on very different timescales. On the one hand, the orbital evolution spreads over millions of years or more, demanding large calculation time steps, usually of the order of years (for instance, our sample of model JFCs has a fixed output frequency of 1 per 100 yr). On the other hand, typical timescales for processes at the origin of cometary activity can vary from a few hours or days (e.g., sublimation of volatile species; De Sanctis et al. 2015), to months or years (e.g., crystallization of amorphous water ice in the giant-planet region; Guilbert-Lepoutre 2012). Solving the time-dependent equations of heat transfer and gas flow in a porous medium, while accounting for multiple phase transitions (e.g., crystallization, sublimation), during an orbital evolution spanning on million years, leads to a prohibitive amount of calculation time. As a consequence, a number of simplifying assumptions are adopted in order to keep the problem at hand tractable.

### 2.1. Thermal Evolution Model

First, a driving assumption comes from the 1 per 100 yr dynamical output frequency. For such a long time step, latitudinal effects at the surface, due to the shape of a comet nucleus, its rotation (diurnal variations), or seasonal variations, cannot be resolved. Because they need to be averaged out during this dynamical time step, there is no critical need to use a slow-rotator approximation for the surface (Huebner et al. 2006) or a 2D/3D thermal evolution model for the interior. Besides, in our sample, the vast majority of clones spend a substantial fraction of their lifetime at large heliocentric distances, on transition orbits between Saturn and Neptune (e.g., Figure 3), where diurnal and seasonal temperature variations are less significant. Therefore, we can consider that a 1D approximation is adequate for this study. We thus use a 1D version of the 3D thermal evolution model described by Guilbert-Lepoutre et al. (2011) to solve the heat diffusion equation:

$$ \rho_{\text{bulk}} c \frac{\partial T}{\partial t} + \text{div}(-\kappa \overrightarrow{\text{grad}} T) = S, \quad (1) $$

where $\rho_{\text{bulk}}$ (kg m⁻³) is the object's bulk density, c (J kg⁻¹ K⁻¹) the material's heat capacity, T (K) the temperature, $\kappa$ (W m⁻¹ K⁻¹) the material's effective thermal conductivity, and S the heat sources and sinks.

Second, we apply simplifying assumptions to the thermophysical properties of cometary material. The effective thermal conductivity of cometary material can be written as

$$ \kappa = h \phi \kappa_{\text{solid}} + \kappa_{\text{rad}}, \quad (2) $$

where h and φ are reduction factors for which various expressions exist (e.g., Shoshany et al. 2002; Gundlach & Blum 2012; Ferrari & Lucas 2016; see Figure 1 for comparison). The factor h is a dimensionless quantity, known as the Hertz correction factor. It

<figure><img src="image_4.png" /><figcaption>Figure 1. Values of the thermal conductivity $\kappa = h \phi \kappa_{\text{solid}}$, with $\kappa_{\text{solid}} = 4.2$ (W m⁻¹ K⁻¹), and correction factors computed with various formulae available in the literature.</figcaption></figure>

describes the reduction of the effective cross section of the grains in a porous material (Gundlach & Blum 2012) and may theoretically vary from $10^{-4}$ to 1, although values of the order of $10^{-2}$ are considered the most plausible and commonly used (see Huebner et al. 2006 for a review). The factor φ is a correction factor applied to account for the effect of the porous structure of cometary material: Here we use Russell's correction factor (Russell 1935), calculated as

$$ \phi = \frac{\psi^{2/3}f + (1 - \psi^{2/3})}{\psi - \psi^{2/3} + 1 - \psi^{2/3}(\psi^{1/3} - 1)f}, \quad (3) $$

where ψ is the porosity and f is the ratio between the solid conductivity and the radiative conductivity $\kappa_{\text{rad}}$, which accounts for the transfer of heat through radiation in the pores:

$$ \kappa_{\text{rad}} = 4r_p \varepsilon \sigma T^3, \quad (4) $$

with $r_p$ (m) the average pore radius, usually set to be the same size as the grains of the medium (Huebner et al. 2006; so ~1 µm in our case), ε the material's emissivity, and σ the Stefan-Boltzmann constant.

In Figure 1, we can see that the various expressions for the thermal conductivity found in the literature produce significantly different values, varying by several orders of magnitude, depending on the porosity or value for the Hertz factor. We note that heat transport by radiation inside the pores becomes significant only for large values of the pore size and at temperatures higher than 100–150 K (higher than 200 K for the largest values of thermal conductivity). We thus ignore $\kappa_{\text{rad}}$ in Equation (2) and use Russell's equation, further corrected by the Hertz factor. Three values—$10^{-2}$, $10^{-3}$, and $10^{-4}$— are considered for this factor in order to assess the influence of thermal conductivity on the simulation outcomes.

Finally, the most crucial simplification is to ignore phase transitions. In Equation (1), this translates to $S = 0$. This
<table><caption>**Table 1**

Thermal Model's Physical Parameters</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>Radius</td><td><i>R</i></td><td>5000</td><td>m</td><td></td></tr><tr><td>Initial temperature</td><td><i>T</i></td><td>10</td><td>K</td><td></td></tr><tr><td>Porosity</td><td><i>ψ</i></td><td>0.8</td><td></td><td>Kofman et al. (2015)</td></tr><tr><td>Bulk density</td><td><i>ρ</i><sub>bulk</sub></td><td>525</td><td>kg m<sup>-3</sup></td><td></td></tr><tr><td>Dust density</td><td><i>ρ</i><sub>dust</sub></td><td>3500</td><td>kg m<sup>-3</sup></td><td>Huebner et al. (2006)</td></tr><tr><td>Heat capacity</td><td><i>c</i><sub>dust</sub></td><td>1000</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>Kömle et al. (2017)</td></tr><tr><td>Effective conductivity</td><td><i>κ</i></td><td>4.2</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>Ellsworth & Schubert (1983)</td></tr><tr><td>Mean pore radius</td><td><i>r</i><sub>p</sub></td><td>10<sup>-6</sup></td><td>m</td><td>Huebner et al. (2006)</td></tr><tr><td>Bond albedo</td><td><i>A</i></td><td>0.04</td><td></td><td>Huebner et al. (2006)</td></tr><tr><td>Emissivity</td><td><i>ε</i></td><td>0.9</td><td></td><td>Huebner et al. (2006)</td></tr><tr><td>Hertz factor</td><td><i>h</i></td><td>10<sup>-2</sup>–10<sup>-4</sup></td><td></td><td>Huebner et al. (2006)</td></tr></tbody></table>

means that our model does not faithfully describe the thermal processing during active stages of the orbital evolution when phase transitions are at the origin of cometary activity. It also entails that no gas diffusion is accounted for, such that the mass is conserved, allowing the use of composition-independent thermophysical properties for the cometary material, in particular for the heat capacity and the thermal conductivity.

Each simulated JFC is thus modeled as a highly porous sphere with a radius of 5 km. Thermophysical parameters are chosen to be averages in the published literature (see Table 1). The initial temperature is set at 10 K, allowing for any type of temperature increase or decrease related to the orbital evolution, and significantly lower than the sublimation temperatures of hypervolatiles. This aspect is important for examining the relative depths where conditions for the loss of such free condensed species are present.

The boundary condition at the surface for Equation (1) is

$$ (1 - A) \frac{L_{\odot}}{4\pi d_{H}^{2}} = \varepsilon \sigma T^{4} + \kappa \frac{\partial T}{\partial z} \quad (5) $$

with the nucleus' insolation given as a function of *A* the Bond's albedo, *L*<sub>☉</sub> the solar constant and *d*<sub>H</sub> (au) the heliocentric distance; the thermal emission given as a function of *ε* the emissivity, *σ* the Stefan–Boltzmann constant and *T* (K) the temperature; and the heat flux toward the interior given as a function of the surface's thermal conductivity *κ* (W m<sup>-1</sup> K<sup>-1</sup>). We assume that the incident solar energy is uniformly distributed over the surface of the nucleus, an approximation known as “fast rotator,” providing a spherical average of the energy received by the nucleus (Huebner et al. 2006).

## 2.2. Model JFCs Sample

Our sample of model JFCs was produced from a simulation by Nesvorný et al. (2017) of the long-term dynamical evolution of outer solar system bodies. It includes 276 model JFCs chosen from a run that successfully reproduced the orbital distribution of observed JFCs. This run was part of a simulation performed with moderate timescales describing the implantation of the Kuiper Belt and scattered disk, assuming a two-stage migration for Neptune: 10 Myr for the first stage, before the instability, and 30 Myr for the second stage after the instability (see Section 3 in Nesvorný et al. (2017), for a detailed description). The sample of model JFCs was produced during the last segment of the integration, from *t* ≈ 3.5 Gyr

until the current epoch. Just like actual JFCs, the simulated ones start in their source reservoir, mainly the scattered disk, and gradually evolve inward during the simulation. All selected model JFCs were active and observable, having at some point in their lifetime perihelion distances within 2.5 au. For each model JFC, the orbital evolution starts to be recorded from the first time its heliocentric distance is within 30 au and continues until its ejection from the solar system. The dynamical lifetimes of our model JFCs range from a few thousand years (minimum value ~0.1 Myr) to almost a billion years (maximum value ~908 Myr), with a median of ~54 Myr. The evolution of orbital elements (heliocentric distance, semimajor axis, eccentricity, and inclination) is recorded every 100 yr.

Figure 2 presents an example of the orbital evolution for a model JFC, from the moment it is first detected within Neptune's orbit, until its ejection. This model JFC has a relatively long lifetime of ~19 Myr. Figure 2 demonstrates the hand-off process through which a comet is scattered inwards from one planet to the planet interior to it until it reaches Jupiter-crossing orbits (Levison & Duncan 1997). At this point the comet's semimajor axis is reduced such that it passes within 2.5 au and becomes an observable JFC (Brasser & Wang 2015; Roberts & Muñoz-Gutiérrez 2021). It is clear from the density of points in Figure 2 that the model JFC spends most of its lifetime (~16 Myr) beyond 10 au, with orbits between those of Neptune and Saturn. During this time, it would be classified as a Centaur, defined by Jewitt (2009) as an object with *a*<sub>J</sub> < *a*, *q* < *a*<sub>N</sub>. Once under the dynamical control of Jupiter, the ejection phase begins, and the model JFC is scattered outwards and quickly leaves the solar system (~3 Myr).

Figure 3 shows where in the orbital parameter space our 276 model JFCs statistically spend their time. As was the case for the model JFC from Figure 2, most model JFCs tend to spend the bulk of their lifetimes, once within 30 au, on transient orbits between those of Saturn and Neptune in the Centaur area. Of course, there is a diversity of outcomes governed by the stochastic nature of close gravitational encounters with the giant planets. JFCs scattered quickly inward by Neptune, Uranus, and Saturn tend to have shorter dynamical lifetimes. Indeed, every model JFC's dynamical lifetime is short once it becomes a true JFC with *q* < 2.5 au (Di Sisto et al. 2009; Nesvorný et al. 2017). From the distribution in Figure 3, we note that there is no clear limit in eccentricity values, although values substantially cluster between 0.1 and 0.5.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example of a model JFC’s orbital evolution from the first time it crosses Neptune’s orbit, until its ejection from the solar system after approximately 19 Myr. The orbits of the giant planets are given by the black dots. Black curves show the locus of orbits with perihelia corresponding to the distance of each giant planet. The color code provides the time evolution.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Density histogram of the orbits for the entire sample in the semimajor axis–eccentricity plane. Orbits of the giant planets are highlighted with white and black dots. The color code represents the number of orbits per bin width.</figcaption></figure>

## 2.3. Coupling the Thermal and Dynamical Evolution

In practice, the thermal and dynamical evolution of each model JFC are coupled via the boundary condition of the thermal evolution model (Equation (5)), where the heliocentric

distance is required to compute the energy balance at the surface. Figure 4 shows an example of the evolution for a model JFC’s fully recorded orbital evolution. We note that the orbits explored span a wide range of orbital periods, sometimes
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of orbital parameters for a model JFC. Upper panel: semimajor axis, perihelion and aphelion distances, and equivalent semimajor axis. Middle panel: eccentricity and inclination of the model JFC. Lower panel: orbital period of the model JFC.</figcaption></figure>

larger than the dynamical time step (100 yr), sometimes much smaller (Figure 4, bottom panel). Hence, only a fraction of the orbit is explored (or in the case of orbits with small orbital periods, several orbits are explored during the dynamical time step, plus a fraction): From a thermal point of view, large differences can arise whenever this fraction is close to perihelion or aphelion. To alleviate this problem, and in order to treat all orbits in the same fashion, we use an energy-averaged orbital distance. Prialnik & Rosenberg (2009) showed that for simulations on long timescales, eccentric orbits can be modeled as circular ones receiving the same total energy over an orbital period. The equivalent semimajor axis $a_{eq}$ is simply $a_{eq} = a(1 - e^2)$ with *a* the semimajor axis and *e* the orbital eccentricity. The evolution of $a_{eq}$ is shown in red in Figure 4, upper panel, between the semimajor axis and perihelion. In our thermal simulations, we recalculate $a_{eq}$ every 100 yr, when orbital elements change. For the duration of the dynamical output, this equivalent circular orbit is subsampled in order to perform thermal calculations on a smaller, more appropriate thermal time step.

## 3. Results

### 3.1. Processing of Individual Model JFCs

Figure 5 shows the evolution of the temperature distribution for the same model JFC as in Figure 4. It is a relatively short-lived one within our sample: Quickly scattered inwards by Neptune, such that it orbits close to Saturn almost from the beginning of the dynamical run (at ~3500 Myr), it is ejected from the solar system after ~1.45 Myr. We note that the evolution of the temperature distributions depends on the

thermal conductivity, depicted here via the Hertz correction factor (Figure 5). Using the 80 K isotherm as a rough indicator for the sublimation of moderately volatile species (such as CO$_2$; white dashed line in Figure 5), we observe that the lower the Hertz factor, the lesser the extent of internal heating: The 80 K isotherm is located at ~150 m for $h = 10^{-2}$, at ~50 m for $h = 10^{-3}$, and at ~10 m for $h = 10^{-4}$. This is a natural consequence of heat conduction, as the skin depth $\delta$—a rough estimate of how deep a surface temperature change can propagate below the surface—is directly dependent on the value of the thermal conductivity ($\delta = \sqrt{\frac{P\kappa}{\pi\rho_{\text{bulk}}c}}$, with P a reference period of time).

Figure 5 also shows how the chaotic orbital evolution of a comet nucleus is imprinted on its thermal evolution. This particular model JFC enters (and then leaves) regions of thermal significance in the inner solar system (e.g., where critical phase transitions such as water-ice sublimation can be triggered) more than once during its lifetime. Three main approaches can be identified in this case, associated with individual heating episodes (most clearly visible in panel (b) of Figure 5: two at $t \sim 3502.85$ Myr and one at $t \sim 3503.3$ Myr). This trend of multiple heating episodes, observed for all model JFCs in our sample, with a varying number of approaches and varying degrees of heating intensity (depending on both their duration and the associated equivalent semimajor axis value), suggests that actual comet nuclei observed in the inner solar system may have been exposed to substantial heating in their past, allowing for thermally induced alteration processes to occur in deep layers below their surface.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Temperature distributions in a subsurface layer of 200 m for the model JFC of Figure 4, and for three different values of the Hertz factor. Panels: (a) Evolution of the equivalent semimajor axis (see Section 2.2), (b) temperature profile evolution over time for *h* = 10⁻², (c) temperature profile evolution over time for *h* = 10⁻³, and (d) temperature profile evolution over time for *h* = 10⁻⁴.</figcaption></figure>

## 3.2. Processing in the Sample of Model JFCs

To constrain statistically the thermal processing of the model JFCs in our sample, we track the depths of three isotherms, at two specific moments of their dynamical lifetime: (a) the moment a model JFC transitions to a JFC orbit, allowing us to examine the initial conditions for JFCs, and (b) the last moment of their JFC phase, in order to study the degree of total processing.

For each phase transition, a characteristic timescale can be computed (see Prialnik et al. 2004 for sublimation timescales), such as $\tau_{cr} = 9.54 \times 10^{-14} \exp\left(\frac{5370}{T}\right)$ for the crystallization of amorphous water ice (Schmitt et al. 1989). In our case, we assume a characteristic timescale of the order of 100 yr (the dynamical time step), which provides a corresponding temperature associated with these common alteration processes:

1. 25 K, representative of sublimation temperatures of hypervolatile species (such as CO);

2. 80 K, representative of sublimation temperatures for moderately volatile species (such as CO₂);

3. 110 K, representative of the amorphous to the crystalline water-ice phase transition, a relevant process as indirect observational data suggest the presence of amorphous ice both in comets (Meech et al. 2009) and Centaurs (Jewitt 2009; Guilbert-Lepoutre 2012).

We adopt the Sarid et al. (2019) classification, where JFCs are defined as objects with $q < 5.2$ au and $Q < 7$ au. We avoid classifications using the Tisserand parameter with respect to Jupiter (Levison 1996): Although more accurate, they ignore the presence of the other planets. In our sample, some model JFCs (~11%) do not satisfy the Sarid et al. (2019) criterion for JFCs. For those, we consider that the transition to a JFC occurs the first time at $q < 2.5$ au.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Location of the 110 K isotherm for the entire sample at $t_0$, the moment of transition to JFC (blue solid line), and $t_f$, the final moment of the JFC phase (orange dashed line) for the three values of the Hertz factor considered: $h = 10^{-2}$ (panel a), $h = 10^{-3}$ (panel b), and $h = 10^{-4}$ (panel c).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Location of the 80 K isotherm for the entire sample at $t_0$, the moment of transition to JFC (blue solid line), and $t_f$, the final moment of JFC phase (orange dashed line) for the three values of the Hertz factor considered: $h = 10^{-2}$ (panel a), $h = 10^{-3}$ (panel b), and $h = 10^{-4}$ (panel c).</figcaption></figure>

For the sake of illustration, we use the model JFC presented in Figure 5. Following the Sarid et al. (2019) criterion, it transitions to a JFC orbit for the first time at $t=3502.8687$ Myr, having previously spent ~0.88 Myr on transient orbits in the giant-planet region. At the moment of this first transition, the three isotherms (25 K, 80 K, 110 K) are located at ~810 m, ~60 m, and ~5 m, respectively, below the surface for $h = 10^{-2}$, and at ~60 m, ~3 m, and ~60 cm, respectively, for $h = 10^{-4}$. At the final moment of its JFC phase, having experienced intense processing, for $h = 10^{-2}$, the 25 K isotherm is located at the same depth (~810 m), whereas the 80 K isotherm advanced at ~160 m and the 110 K isotherm at ~60 m. For $h = 10^{-4}$, the three isotherms are closer to the surface (at ~60 m, ~9 m, and ~4 m, respectively). We emphasize that because no phase transition is actually computed, there is no heating delay associated with the sublimation of volatile species nor any recondensation. Indeed, when volatile species sublimate, their partial pressure peaks at

the phase transition front. Molecules then follow pressure gradients and recondense in any volume that sustains the appropriate temperature and pressure conditions, either toward the surface or deeper in the interior. We track the temperature evolution for timescales so long that both sublimation and gas diffusion timescales are shorter than the lifetime of our model JFCs. As such, the results presented here and below are significant with respect to the primordial composition, and whether it can be maintained in near-surface layers. We thus assume that volatile species, initially present as pure condensate, can sublimate if their sublimation temperature is reached: free condensed CO would sublimate down to the level of the 25 K isotherm, and free condensed CO₂ down to the level of the 80 K isotherm. Whether molecules escape the nucleus or recondense in the interior, where conditions allow it, the net result is an alteration of the primordial composition. When these layers subsequently contribute to the cometary activity that we observe, they are no longer pristine, thus they
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Location of the 25 K isotherm for the entire sample at $t_0$, the moment of transition to JFC (blue solid line), and $t_f$, the final moment of the JFC phase (orange dashed line) for the three values of the Hertz factor considered: $h = 10^{-2}$ (panel a), $h = 10^{-3}$ (panel b), and $h = 10^{-4}$ (panel c).</figcaption></figure>

do not reflect the primitive material incorporated in comet nuclei. For simplicity, we refer to the “loss” of respective volatile species hereafter: The reader should keep in mind that it is the primordial inventory of volatile species, which is lost.

With the above assumptions, the above model JFC would have “lost” these volatile species from a substantial subsurface layer before its transition to JFC. During the intense processing of the JFC phase, we note that it is the moderately volatile species and the crystallization of amorphous water ice that are mostly considered, as the fate of its hypervolatile content has already been determined by the time spent in transient orbits at the Centaur area. We present below the potential processing for the whole sample.

### 3.2.1. Crystallization of Amorphous Water Ice

We assume here that the crystallization of amorphous water ice, if at all present inside JFCs, would follow the 110 K isotherm. Figure 6 shows the depths of this isotherm, as a function of the Hertz factor, for the initial and final moments of the model JFCs. Before transitioning to a JFC orbit for the first time, the three different values of the conductivity yield similar depth distributions: For nearly all model JFCs in the sample, the 110 K isotherm remains located within 40 m below the surface (up to 80 m for $h = 10^{-2}$). As can be expected from the conduction of heat with time, the isotherm propagates below the surface during the JFC phase (for all values of the effective thermal conductivity), remaining, however, quite close to the surface: For the bulk of the model JFCs, the isotherm is located in the top ~80 m below the surface. In the case of $h = 10^{-2}$, we can find it in deeper layers, though rarely below ~400 m. These depths are compatible with the results from Guilbert-Lepoutre (2012), who studied the survival of amorphous water ice in Centaurs, using fixed orbits in the giant-planet region for 10 Myr. Our results suggest that the survival of amorphous water ice is definitely conceivable in the interior of comet nuclei and do not challenge the hypothesis that this phase transition to crystalline water ice could be a possible source for activity among Centaurs or possible outbursts (e.g., Wierzchos & Womack 2020).

### 3.2.2. Volatile Content Following the 80K Isotherm

The distributions of depths reached by the 80 K isotherm are quite similar to those for the 110 K isotherm (see Figure 7), although because the temperature of interest is smaller, it is generally found deeper below the surface. The minimum processing of model JFCs in our sample involves the top ~100 m for most objects and all values of thermal conductivity. For $h = 10^{-2}$, we note, however, that almost all model JFCs are liable to lose their moderately volatile content in the top 40 m below the surface, as the 80 K isotherm has progressed below that limit even before the model JFCs transition to JFC orbits: 69.5% of the model JFCs in the sample have the isotherm located between 40 and 80 m, 28% have it below 80 m. Only 2.5% (seven model JFCs) have the 80 K isotherm above 40 m. A similar trend is observed for $h = 10^{-3}$, albeit less pronounced: 62.3% of model JFCs have the 80 K isotherm between 40 and 80 m in this case. For $h = 10^{-4}$, the majority of the model JFCs (79%) have this isotherm located in the first 40 m (average depth of 17.2 m).

At the end of their JFC phase, the isotherm is found ~80 m below the surface on average for $h = 10^{-2}$, ~70 m for $h = 10^{-3}$, and ~9 m for $h = 10^{-4}$. This suggests that the dominant heating phase related to this particular isotherm occurs during their JFC phase, although the processing during the Centaur phase is not negligible at all, even for the least conductive scenario. Accounting for the sublimation of the corresponding volatile species on long timescales will therefore be of particular importance for future works.

### 3.2.3. Hypervolatile Content Traced by the 25K Isotherm

The distribution of depths reached by the 25 K isotherm for the different values of the thermal conductivity is given in Figure 8. It is apparent that with respect to this particular isotherm, the thermal conductivity (through the choice of Hertz factor in our study) is critical to assess the corresponding processing, while the time at which the processing is examined plays no significant role. Indeed, the distributions obtained at the first time of transition and at the end of the JFC phase are almost identical (especially in the case $h = 10^{-2}$). For $h = 10^{-2}$
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Maximum depths of the 25 K isotherm as a function of the lifetimes of the 276 model JFCs of our sample for three different values of the Hertz factor.</figcaption></figure>

(and $10^{-3}$ to some degree), most model JFCs have already heated up down to their cores before their first transition to a JFC orbit (77% and 33.3%, respectively). For the least conductive scenario ($h = 10^{-4}$), only 27.9% of model JFCs retain the 25 K isotherm in the top 250 m layer, the average depth for the overall sample being ~650 m. At the final moment of their JFC lifetime, we see that less than 10% of model JFCs are able to maintain the 25 K isotherm in the top 1 km (only one model JFC keeps the isotherm in the top 250 m) for $h = 10^{-2}$. For $h = 10^{-3}$, only 30% of model JFCs keep the isotherm within 1 km below the surface (~8% keep it within 250 m), whereas for $h = 10^{-4}$, 27% of model JFCs have the 25 K isotherm within 250 m, and 50% within 500 m.

Assuming that the 25 K isotherm is indicative of the loss of CO, for example, all model JFCs in our sample (except one) would have lost free condensed CO in the top 250 m subsurface layer before their first arrival on JFC orbits. We could further argue that, given the extremely low sublimation temperature of hypervolatiles, these species would be lost as pure condensate even before the model JFCs enter the giant-planet region because isothermal surface temperatures in the Kuiper Belt range between 30 and 50 K. Because the parameter controlling the outcomes of thermal simulations with respect to the 25 K isotherm is the thermal conductivity, the time spent in transient orbits is also critical. We illustrate this effect in Figure 9: Model JFCs spending long periods of time on transient orbits (i.e., model JFCs with longer lifetimes) are more susceptible to being heated up at greater depths. For $h = 10^{-2}$, model JFCs with lifetimes over ~8 Myr are prone to losing all free condensed hypervolatiles down to their core. For $h = 10^{-3}$, model JFCs with lifetimes above ~90 Myr would experience the same fate. In contrast, no model JFC is found to lose free condensed hypervolatiles below 2500 m for $h = 10^{-4}$. The exact depth will need to be confirmed by including the relevant

phase transition, with appropriate diurnal and seasonal approximations, in further modeling work.

### 3.3. Relation to the Observed JFC Population

We use a random number generator to select a moment during the JFC phase of all model JFCs in order to examine the thermal processing of typically observed JFCs. The depth distributions for the different temperatures and Hertz factors considered are presented in Figure 10. We find that the depth distributions for the 25 K isotherm are almost identical to the distributions during the transition to JFCs and the JFC end states, with only slight differences of the order of 2%–3% for model JFCs having the isotherm close to the surface. This confirms our previous argument that for pure condensed hypervolatiles, the key factor is the time spent prior to the JFC phase in the Centaur area or before and the thermal characteristics of the cometary material controlling the amount of heat transferred toward the interior of comet nuclei.

As expected, the 80 K isotherm is located slightly deeper below the surface (for all thermal conductivities) when compared to the moment of transition and slightly shallower when compared to the final states of JFCs. For instance, in the case of $h = 10^{-2}$, although in 38% of the model JFCs the isotherm is retained between 40 and 80 m, in the remaining 62% it is located significantly deeper, as far as ~750 m, whereas no model JFC has the isotherm in the first 40 m below the surface. Similar behavior is observed for $h = 10^{-3}$ and $h = 10^{-4}$, as witnessed by the increase in the average depths presented in Table 2. Similar results are drawn for the 110 K isotherm: For a randomly selected moment of the JFC phase, it is located between the initial and the final positions shown in Figure 6. Taking for example the case of $h = 10^{-2}$, we see that at a random moment, the fraction of model JFCs having this isotherm beyond the first 80 m below the surface is ~30%, whereas in the moment of transition, this fraction was
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Temperature distributions for the entire sample at a random moment of their JFC phase for the different values of the Hertz factor considered: *h* = 10⁻² (blue solid line), *h* = 10⁻³ (orange dashed line), and *h* = 10⁻⁴: (a) 25 K and (b) 80 K (110 K).</figcaption></figure>

less than 10% and at the final moment of the JFC phase mounts to ~40%.

# 4. Discussion

## 4.1. Composition Effects and Latent Heat Transfer

To simplify our study, we have assumed that our model JFCs are made of a porous structure of refractories, with no ices included so as to avoid phase transitions. Despite the uncertainties regarding their structure and composition, actual comets clearly contain a number of icy species. Ignoring their presence removes the energy sources and sinks related to their phase transitions, such as the energy released during the crystallization of amorphous water ice (Schmitt et al. 1989) and the energy consumption during ice sublimation (Prialnik et al. 2004). In fact, when sublimation and recondensation of volatile compounds do occur inside comet nuclei, heat transport via the vapor phase is very important, sometimes more effective than heat conduction by the solid matrix (see Huebner et al. 2006 for a review). As such, during the active phase, the dominant heat transport mechanism is the transfer of latent heat: This strongly influences the temperature profile within actual comet nuclei and is not accounted for in our model JFCs. The direct consequence is that while an icy compound is sublimating, the input energy is used for the phase transition so that the internal temperature will stop increasing. This would limit the penetration of successive heat waves observed in Figure 5, as the heat would be consumed to sublimate ices, rather than propagating inward. In this regard, the depths achieved in our study could be considered as upper limits to the actual processing of comet nuclei.

Ignoring phase transitions also makes it impossible to account for surface erosion, which is associated with water-ice sublimation. Estimates from the literature indicate that the erosion rate can vary from tenths of centimeters up to approximately 2 m per year, depending on thermal characteristics, the heliocentric distance, the inclination of the spin axis, and the cometocentric latitude examined. Huebner et al. (2006) estimate an erosion rate varying from 10 cm yr⁻¹ to 2.1 m yr⁻¹ for 46P/Wirtanen, while estimates for 67P/Churyumov-Gerasimenko give an average erosion rate ranging from 0.67 to 2.9 m, depending on the thermophysical

<table><caption>**Table 2**
Average Depth (m) of Three Isotherms for Three Values of the Hertz Factor, at Distinct Times of Our Sample of Model JFCs' Orbital Evolution</caption><thead><tr><th>T (K)</th><th>h = 10<sup>-2</sup></th><th>h = 10<sup>-3</sup></th><th>h = 10<sup>-4</sup></th></tr></thead><tbody><tr><td colspan="4">First Transition to a JFC Orbit</td></tr><tr><td>25</td><td>4142.2 ± 1529.9</td><td>2561.5 ± 1894.0</td><td>651.7 ± 615.2</td></tr><tr><td>80</td><td>125.7 ± 92.3</td><td>50.1 ± 36.8</td><td>17.2 ± 23.0</td></tr><tr><td>110</td><td>27.1 ± 37.5</td><td>10.0 ± 19.4</td><td>2.1 ± 5.4</td></tr><tr><td colspan="4">Random Time in the JFC Phase</td></tr><tr><td>25</td><td>4156.3 ± 1502.9</td><td>2566.60 ± 1889.8</td><td>654.1 ± 613.5</td></tr><tr><td>80</td><td>209.2 ± 130.1</td><td>81.0 ± 38.4</td><td>31.2 ± 26.9</td></tr><tr><td>110</td><td>89.9 ± 39.9</td><td>32.3 ± 32.8</td><td>10.1 ± 18.0</td></tr><tr><td colspan="4">End of the JFC Phase</td></tr><tr><td>25</td><td>4181.7 ± 1474.1</td><td>2575.1 ± 1882.8</td><td>656.2 ± 611.8</td></tr><tr><td>80</td><td>278.3 ± 177.6</td><td>99.7 ± 51.1</td><td>40.9 ± 28.1</td></tr><tr><td>110</td><td>107.9 ± 103.4</td><td>40.2 ± 39.0</td><td>13.4 ± 21.9</td></tr></tbody></table>

parameters considered for the surface (Keller et al. 2015). Accounting for the cumulative effect of erosion over long timescales becomes important during the active phase of JFCs, as a significant amount of processed surface material can be removed, revealing unprocessed layers (or bringing them closer to the surface) and changing the internal stratigraphy. A more sophisticated model, including carefully developed approximations for diurnal and seasonal variations of the temperature and activity, is clearly required for quantitative calculations of the actual internal stratification or the cumulative erosion with time.

Finally, we have seen how the thermal conductivity is instrumental in constraining the extent of subsurface thermally induced processing. Arguments summarized in Huebner et al. (2006) are in favor of a Hertz factor of the order of 10⁻² (based on laboratory experiments, observations, and theoretical considerations), which in our simulations corresponds to the maximum processing considered. We have also ignored heat transport through radiation within the porous structure. For the lowest values of the Hertz factor, accounting for $κ_{rad}$ would
actually increase the effective thermal conductivity, even at low temperatures ($\kappa_{solid}$ and $\kappa_{rad}$ are of the same order of magnitude around 40–50 K for the lowest values of the Hertz factor; Huebner et al. 2006). This makes our results for $h = 10^{-4}$ unrealistic and favors evolution outcomes given for $h = 10^{-2}$ or $h = 10^{-3}$, albeit still limited by the lack of icy components in the material.

## 4.2. Limitations of the Orbital Approximation

To couple the thermal model to dynamical simulations, we effectively placed each model JFC on a circular orbit receiving the same total energy as its true, eccentric orbit (see Section 2.2). In practice, this averaging means that a model JFC receives a constant amount of solar energy for every thermal time step, regardless of its true orbit. The input flux and, therefore, the calculated surface temperature (Equation (5)) are both underestimated when actual heliocentric distances are smaller than the equivalent semimajor axis (typically close to the perihelion where the most intensive heating takes place) and overestimated when they are larger, especially close to the aphelion. In this regard, the effects of cumulative perihelion passages are always underestimated in our study. On the other hand, one could argue the model JFCs are not allowed to cool down when they are close to aphelion. After further analysis using real orbits, we find that this approximation is valid for relatively low eccentricities (e.g., $e < 0.5$), where deviations from circular orbits are small, as well as for orbits with a large semimajor axis because temperature variations across the orbit are less significant. This is the case for the majority of the orbits in our sample (see Figure 3). If improvements are required in this study, they need to address the key points identified above, i.e., accounting for phase transition or better constraints on the thermophysical parameters, rather than on the orbital approximations.

## 4.3. On the Activity of Centaurs and JFCs

With the aforementioned limits in mind, our results have several inferences regarding the activity of Centaurs and JFCs. Before our model JFCs go through their first transition to a JFC orbit, we consider them as Centaurs. Recent studies have tried to link the activity of Centaurs to their orbital evolution. Indeed, Guilbert-Lepoutre (2012) found that crystallization of amorphous water ice can be a source of activity, efficient at heliocentric distances up to 10–12 au, sustained though for a limited time (typically hundreds to thousands of years). This implies that active Centaurs should have suffered a recent orbital change. Davidsson (2021) nuanced this result by showing that crystallization would be efficient only up to 8–10 au, while the sublimation of CO$_2$ would be the only source of activity in the 10–12 au region. The study by Fernández et al. (2018) is compatible with a link between thermal and orbital evolution as the origin of active Centaurs. They find that active Centaurs are more prone to drastic drops in their perihelion distances than inactive Centaurs, with timescales of the order of $10^2–10^3$ yr. Finally, Cabral et al. (2019), Li et al. (2020), and Lilly et al. (2021) found no activity among various Centaur detections: The corresponding objects are dynamically stable on long timescales.

Our study shows that the 25 K isotherm reaches deep layers below the surface of most model comets early in the Centaur phase. Indeed, heating subsurface layers above this 25K temperature starts at large heliocentric distances and is relevant

throughout the entire orbital evolution. Both the sublimation of moderately volatiles (such as CO$_2$) and the crystallization of amorphous water ice would be triggered during this Centaur phase, each time their perihelion distance drops and the heat wave reaches an internal layer that was not previously depleted during its past orbital evolution. Our sample suggests that on average, our objects could be active for a $10^3$ period during their $10^6–10^7$ yr life as Centaurs. This period is not continuous: Following chaotic evolution in the giant-planet region, we see that our model JFCs are active for several shorter periods of time, typically ~10 times. We emphasize nonetheless that our results cannot be directly compared to the observed Centaur population because both the latter and our dynamical sample are biased. On the one hand, there has not been any systematic survey targeting the Centaur population (e.g., Cabral et al. 2019), so the typical ~10% fraction of active Centaurs is meant to evolve. Besides, active Centaurs are more easily observable, typically being closer to the Sun, and displaying cometary activity. On the other hand, all our model Centaurs/JFCs do become JFCs at one point of their lifetime, which is not the case for all actual Centaurs, as some may be directly ejected before reaching Jupiter-crossing orbits.

## 4.4. Considerations on Observed JFCs

Following the thermal processing during the Centaur phase, once our model JFCs transition to actual JFC orbits, their subsurface layers can be already considerably altered. Our results suggest that at a random moment of their JFC phase, regardless of the temperature examined and the conductivity scenario considered, all model JFCs have already undergone sufficient heating for alteration processes to take place at considerable depths. We can thus infer that a typical observed JFC can be significantly altered, with the primordial inventory of volatiles lost due to sublimation, gas diffusion, depletion, and enrichment of internal volumes on long timescales. In particular, the processed layer is substantially larger than the estimated size of the layers involved in producing cometary activity observed for most JFCs (e.g., Huebner et al. 2006). In other words, cometary activity as it is observed today should be produced from layers that have been substantially processed. Our results suggest nonetheless that if these layers may be able to retain some amorphous water ice on the one hand, they should on the other hand have lost all pure primitive condensed hypervolatiles and a fraction of their primitive moderately volatile content. Even if we were able to observe an object transitioning to a JFC orbit—while actually knowing it is the first time the object experiences this transition— our results imply that we would be observing an altered body, in particular one with severe modifications of its hypervolatile content. This highlights the importance of accounting for the past history of each comet in order to better constrain its present thermal, physical, and chemical state.

## 4.5. Implications for Future Cryogenic Sample Return Missions

Our study, despite being a first-order approximation of the thermal processing of JFCs, can have direct implications for future cryogenic sample return missions. We observed that the expected loss of hypervolatile ices should be severe from a substantial subsurface area expanding several thousand meters. The same can be said for moderately volatile species, although the subsurface volume concerned is greatly reduced (see also
Table 2). Sample return missions aiming at depths of 3 m below the surface (e.g., Bockelée-Morvan et al. 2021) would likely not be able to access a primitive layer containing pure hypervolatile or moderately volatile ices. As a consequence, cryogenic temperatures lower than 80 K should not be a critical constraint for such missions.

## 5. Summary

We present the results of a first-order study of the coupled thermal and dynamical evolution of JFCs on their path inward from the Kuiper Belt and scattered disk. We applied a simplified thermal evolution model to a sample of 276 model JFCs, taken from a dynamical simulation that successfully reproduces their observed orbital distribution (Nesvorný et al. 2017). Our simulations show that:

1. Comet nuclei undergo multiple heating episodes, spread randomly across their dynamical lifetimes. These heating episodes are prolonged long-lasting periods, with a large number of orbital changes, each orbit being further characterized by seasonal cycles between their perihelion-aphelion passages. This pattern is observed for all model JFCs during their chaotic transition toward the inner solar system.

2. As a consequence, a substantial subsurface layer is heated, providing the necessary conditions for extensive thermal processing to occur. Processed layers can extend as deep as ~4100 m on average for temperatures allowing the sublimation of hypervolatile species, ~125 m for temperatures permitting the sublimation of moderately volatile species, and ~27 m for temperatures allowing the crystallization of amorphous water ice (considering the most plausible scenario for thermal conductivity). These results have direct implications for the drilling depths of any cryogenic sample return mission, although a more detailed model would be required to confidently conclude limiting depths and temperatures.

3. Despite the limitations of our approach, the fate of hypervolatiles is so extreme that we can infer that all primordial condensed hypervolatiles should be lost from layers that subsequently contribute to any observed cometary activity.

4. For any typical observed JFC, activity is very likely triggered from layers that have been thermally processed and lost their primitive inventory of volatiles. This indicates that JFCs are probably inadequate targets for cryogenic sample return missions. It also points to the necessity of taking into consideration their entire evolutionary history, both thermal and dynamical, when interpreting current observations in a broader context.

We would like to thank the referee for a careful and thorough report on our manuscript. This study is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 802699). We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. S.N.R. is grateful to the CNRS's PNP program (Programme National de Planétologie).

Facilities: PSMN, ENS de Lyon.

## ORCID iDs

Anastasios Gkotsinas https://orcid.org/0000-0002-1611-0381

Aurélie Guilbert-Lepoutre https://orcid.org/0000-0003-2354-0766

Sean N. Raymond https://orcid.org/0000-0001-8974-0758

David Nesvorny https://orcid.org/0000-0002-4547-4301

## References

Bockelée-Morvan, D., Filacchione, G., Altwegg, K., et al. 2021, *ExA*, in press

Brasser, R., & Morbidelli, A. 2013, *Icar*, 225, 40

Brasser, R., & Wang, J. H. 2015, *A&A*, 573, A102

Cabral, N., Guilbert-Lepoutre, A., Fraser, W. C., et al. 2019, *A&A*, 621, A102

Capria, M. T., Coradini, A., De Sanctis, M. C., & Orosei, R. 2000, *AJ*, 119, 3112

Choi, Y.-J., Cohen, M., Merk, R., & Prialnik, D. 2002, *Icar*, 160, 300

Davidsson, B. J. R. 2021, *MNRAS*, 505, 5654

de la Fuente Marcos, C., de la Fuente Marcos, R., Licandro, J., et al. 2021, *A&A*, 649, A85

De Sanctis, M. C., Capaccioni, F., Ciarniello, M., et al. 2015, *Natur*, 525, 500

De Sanctis, M. C., Capria, M. T., & Coradini, A. 2001, *AJ*, 121, 2792

De Sanctis, M. C., Capria, M. T., Coradini, A., & Orosei, R. 2000, *AJ*, 120, 1571

Di Sisto, R. P., Fernández, J. A., & Brunini, A. 2009, *Icar*, 203, 140

Ellsworth, K., & Schubert, G. 1983, *Icar*, 54, 490

Farnham, T. L., Kelley, M. S. P., & Bauer, J. M. 2021, *PSJ*, 2, 236

Fernández, J. A., Helal, M., & Gallardo, T. 2018, *P&SS*, 158, 6

Ferrari, C., & Lucas, A. 2016, *A&A*, 588, A133

Guilbert-Lepoutre, A. 2012, *AJ*, 144, 97

Guilbert-Lepoutre, A., Lasue, J., Federico, C., et al. 2011, *A&A*, 529, A71

Gundlach, B., & Blum, J. 2012, *Icar*, 219, 618

Huebner, W. F., Benkhoff, J., Capria, M.-T., et al. 2006, Heat and Gas Diffusion in Comet Nuclei (Bern: International Space Science Institute)

Hui, M.-T., Farnocchia, D., & Micheli, M. 2019, *AJ*, 157, 162

Hui, M.-T., Jewitt, D., & Clark, D. 2018, *AJ*, 155, 25

Jewitt, D. 2009, *AJ*, 137, 4296

Jewitt, D., Hui, M.-T., Mutchler, M., et al. 2017, *ApJL*, 847, L19

Keller, H. U., Mottola, S., Davidsson, B., et al. 2015, *A&A*, 583, A34

Kofman, W., Herique, A., Barbin, Y., et al. 2015, *Sci*, 349, 2.639

Kömle, N. I., Macher, W., Tiefenbacher, P., et al. 2017, *MNRAS*, 469, S2

Levison, H. F. 1996, in ASP Conf. Ser., 107, Completing the Inventory of the Solar System, ed. T. Rettig & J. M. Hahn (San Francisco, CA: ASP), 173

Levison, H. F., & Duncan, M. J. 1997, *Icar*, 127, 13

Li, J., Jewitt, D., Mutchler, M., Agarwal, J., & Weaver, H. 2020, *AJ*, 159, 209

Lilly, E., Hsieh, H., Bauer, J., et al. 2021, *PSJ*, 2, 155

Lin, H. W., Chen, Y. T., Lacerda, P., et al. 2014, *AJ*, 147, 114

Mazzotta Epifani, E., Dotto, E., Ieva, S., et al. 2018, *A&A*, 620, A93

Mazzotta Epifani, E., Perna, D., Dotto, E., et al. 2017, *A&A*, 597, A59

Meech, K. J., Kleyna, J. T., Hainaut, O., et al. 2017, *ApJL*, 849, L8

Meech, K. J., Pittichová, J., Bar-Nun, A., et al. 2009, *Icar*, 201, 719

Meech, K. J., & Svoren, J. 2004, in Using Cometary Activity to Trace the Physical and Chemical Evolution of Cometary Nuclei, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tuscon, AZ: Univ. Arizona Press), 317

Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, *ApJ*, 845, 27

Prialnik, D., Benkhoff, J., & Podolak, M. 2004, Modeling the Structure and Activity of Comet Nuclei (Tuscon, AZ: Univ. of Arizona Press), 359

Prialnik, D., & Merk, R. 2008, *Icar*, 197, 211

Prialnik, D., & Rosenberg, E. D. 2009, *MNRAS*, 399, L79

Raymond, S. N., Kaib, N. A., Armitage, P. J., & Fortney, J. J. 2020, *ApJL*, 904, L4

Roberts, A. C., & Muñoz-Gutiérrez, M. A. 2021, *Icar*, 358, 114201

Russell, H. 1935, *J. Am. Ceram. Soc.*, 18, 1

Sarid, G., Volk, K., Steckloff, J. K., et al. 2019, *ApJL*, 883, L25

Schmitt, B., Espinasse, S., Grim, R. J. A., Greenberg, J. M., & Klinger, J. 1989, in Physics and Mechanics of Cometary Materials, ed. J. J. Hunt & T. D. Guyenne, Vol. 302 (Paris: European Space Agency), 65

Shoshany, Y., Prialnik, D., & Podolak, M. 2002, *Icar*, 157, 219

Steckloff, J. K., Sarid, G., Volk, K., et al. 2020, *ApJL*, 904, L20

Wierzchos, K., & Womack, M. 2020, *AJ*, 159, 136

Yang, B., Jewitt, D., Zhao, Y., et al. 2021, *ApJL*, 914, L17