# On the Early Thermal Processing of Planetesimals during and after the Giant Planet Instability

Anastasios Gkotsinas, David Nesvorný, Aurélie Guilbert-Lepoutre, Sean N.

Raymond, Nathan Kaib

► To cite this version:

Anastasios Gkotsinas, David Nesvorný, Aurélie Guilbert-Lepoutre, Sean N. Raymond, Nathan Kaib. On the Early Thermal Processing of Planetesimals during and after the Giant Planet Instability. The Planetary Science Journal, 2024, 5 (11), pp.243. 10.3847/PSJ/ad7f4e. hal-05353432

HAL Id: hal-05353432

https://hal.science/hal-05353432v1

Submitted on 7 Nov 2025

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

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

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

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

# On the Early Thermal Processing of Planetesimals during and after the Giant Planet Instability

Anastasios Gkotsinas<fnref n="1" />, David Nesvorný<fnref n="2" />, Aurélie Guilbert-Lepoutre<fnref n="1" />, Sean N. Raymond<fnref n="3" />, and Nathan Kaib<fnref n="4" />

## Abstract

Born as ice-rich planetesimals, cometary nuclei were gravitationally scattered onto their current orbits in the Kuiper Belt and the Oort Cloud during the giant planets' dynamical instability in the early stages of our solar system's history. Here, we model the thermal evolution of planetesimals during and after the giant planet instability. We couple an adapted thermal evolution model to orbital trajectories provided by N-body simulations to account for the planetesimals' orbital evolution, a parameter so far neglected by previous thermal evolution studies. Our simulations demonstrate intense thermal processing in all planetesimal populations, mainly concerning the hypervolatile ice content. Unlike previous predictions, we show that hypervolatile survival was possible in a significant number of planetesimals of the Kuiper Belt and the Oort Cloud. Planetesimals ejected into interstellar space proved to be the most processed, while planetesimals ending in the Oort Cloud were the least processed population. We show that processing differences between populations are a direct consequence of their orbital evolution patterns, and that they provide a natural explanation for the observed variability in the abundance ratios of CO on cometary populations and on the recent observations of long-distance CO-driven activity on inbound long-period comets.

*Unified Astronomy Thesaurus concepts: Long period comets (933); Short period comets (1452); Planetesimals (1259); Comet dynamics (2213); Comet nuclei (2160); Comet volatiles (2162); Comets (280); Interstellar objects (52)*

## 1. Introduction

Comets have captivated the interest of the scientific community for decades due to their presumed pristine nature and their ability to provide valuable information on the conditions of the primordial solar system and its subsequent evolution (e.g., D. Prialnik & M. Podolak 1995; D. Bockelée-Morvan et al. 2015; M. C. De Sanctis et al. 2015; K. Wierzchos & M. Womack 2020).

On a crude outline of a comet's lifetime, we can identify three main evolutionary phases (K. J. Meech & J. Svoren 2004). First, future cometary nuclei formed as ice-rich planetesimals in the outer parts of the Sun's gas-dominated protoplanetary disk (J. B. Simon et al. 2022). Their orbits evolved under the effects of gas drag (I. Adachi et al. 1976) and the growing planets. The giant planets' dynamical instability destabilized the orbits of cometary nuclei, leading them to be strongly gravitationally scattered and, in large part, ejected from the solar system. A fraction of “lucky” planetesimals were trapped in two outer solar system reservoirs: the Kuiper Belt and its extension, the scattered disk (hereafter KB/SD; H. F. Levison et al. 2008), and the Oort Cloud (hereafter OC; e.g., L. Dones et al. 2015; D. Nesvorný 2018; N. A. Kaib & K. Volk 2022). The second phase of the evolution of comets involves a prolonged residence in the reservoirs on relatively stable orbits for timescales of the order of 1 billion yr. Due to the prevailing low temperatures in these faraway reservoirs,

planetesimals are generally considered to have remained unaltered, conserving their primordial physical and chemical characteristics and justifying their characterization as the most primitive objects of our solar system (P. Weissman et al. 2020). In the third phase, these planetesimals return to the inner solar system<fnref n="5" /> and can then be called “comets.” During this phase, comets are usually detected and observed as they undergo extensive thermal processing due to their gradual approach to the Sun. This final phase almost universally culminates with ejection into interstellar space.

The view that comets are “pristine” and largely unaltered has come into question. A number of recent studies have suggested that thermal processing, and the resulting activity and alterations, may have started during the early phases of comets’ evolution, during their prolonged reservoir residence (e.g., J. Haruyama et al. 1993; D. Prialnik & M. Podolak 1995; M. C. De Sanctis et al. 2001; Y.-J. Choi et al. 2002; Q. Kral et al. 2021; A. Parhi & D. Prialnik 2023), or before reaching the inner solar system, in the giant planet area (A. Guilbert-Lepoutre 2012, 2014; D. Jewitt et al. 2021; A. Gkotsinas et al. 2022; N. A. Kaib 2022). This claim is supported by numerous observations of long-distance activity on Centaurs (e.g., D. Jewitt 2009) or inbound long-period comets (LPCs; e.g., K. J. Meech et al. 2017; M.-T. Hui et al. 2019; T. L. Farnham et al. 2021). Some studies (S. N. Raymond et al. 2020; B. J. R. Davidsson 2021; C. M. Lisse et al. 2022) have even suggested processing at the very early stages of comets' lifetimes, within the planetesimal disk and during their implantation into cometary reservoirs. All

<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](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Laboratoire de Géologie de Lyon: Terre, Planétes, Environnement, CNRS, UCBL, ENSL, F-69622, Villeurbanne, France</footnote>

<footnote n="2">²Department of Space Studies, Southwest Research Institute, 1050 Walnut Street, Suite 300, Boulder, CO 80302, USA</footnote>

<footnote n="3">³Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, F-33615 Pessac, France</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, University of Oklahoma, 440 West Brooks Street, Norman, OK 73019, USA</footnote>

<footnote>Received 2024 May 7; revised 2024 September 10; accepted 2024 September 23; published 2024 November 13</footnote>

<footnote n="5">⁵ Here the terms “inner” and “outer” describe the areas within and beyond the orbit of Neptune, respectively.</footnote>
these studies, even if they do not question the primitive status of comets, point to the necessity of reassessing it, in order to better interpret current observations and improve theoretical models and predictions.

In this paper, we reevaluate the primitive nature of comets. We propose a theoretical framework to model the early thermal processing of ice-rich planetesimals that eventually either are ejected from the solar system or return as Jupiter-family comets (JFCs) or LPCs during their implantation into the KB/SD and the OC. Our framework consists of an adapted thermal and ice evolution model coupled to orbital trajectories of implantation or ejection obtained from *N*-body simulations. These allow us to couple the effects of thermal evolution with the comets' orbital evolution, an aspect that has been neglected in most previous studies. By using a large number of particles, we provide a statistically significant first-order assessment of the processing sustained to the three planetesimal populations (KB/SD, OC, and ejected planetesimals) during the first phase of their lifetime, allowing us to make predictions regarding their future evolution in the reservoirs and on their return to the inner solar system.

In Section 2, we provide a brief description of the *N*-body simulations and the planetesimal samples used in this work. In Section 3, we detail the adapted thermal and ice evolution models developed for this study, and we explain how the coupling between these models and the *N*-body simulations is achieved. In Section 4, we present the results of our simulations, and in Section 5, we discuss their implications on the processing of the three studied populations while analyzing the limitations of our approach and their consequences. A brief summary of the main points of this work is provided in Section 6.

## 2. *N*-body Simulations

We used orbital trajectories obtained from two different sets of *N*-body simulations. Both sets simulated the orbital evolution of planetesimals during the giant planet instability (D. Nesvorný 2018) and also included torques from the Galactic tide (J. Heisler & S. Tremaine 1986) in order to enable the capture of planetesimals in the OC. Both sets of simulations invoked the existence of an additional ice giant planet that was ejected during the instability (D. Nesvorný 2011; D. Nesvorný & A. Morbidelli 2012), and both put stringent limits on the simulations that were deemed successful in terms of their matches to the present-day planets' orbits (see, for example, M. S. Clement et al. 2021). The first set of simulations, described in detail in D. Nesvorný et al. (2017) and D. Vokrouhlický et al. (2019), provided two samples of planetesimals that eventually populated the KB/SD and the OC. The second set of simulations, detailed in S. N. Raymond et al. (2020), provided a sample of planetesimals that were ejected out of the solar system. In the following paragraphs, we briefly describe the main considerations and initial conditions of these simulations and give a detailed description of the samples. We refer the reader to the aforementioned publications for a more comprehensive review.

Both sets of *N*-body simulations tracked the orbital evolution of the four giant planets, placed initially in a multiresonant setup (A. Morbidelli et al. 2007), and of an additional ice giant, introduced to increase the probability of reproducing the solar system's current orbital architecture (D. Nesvorný & A. Morbidelli 2012). A number of massless particles ($10^6$ in the first

and $10^3$ in the second simulation) were placed on low-eccentricity ($<0.1$) and low-inclination ($<0.05$) orbits on an outer disk, extending from a heliocentric distance of 24 to 30 au in the first set of simulations and from 21.4 to 30 au in the second set of simulations. The initial internal time step was set at 0.5 yr for the first simulation and at 24 days for the second one. The influence of terrestrial planets was not considered in any simulation, nor were the effects of nongravitational forces.

The simulations started from the moment the giant planets' instability was triggered. The first was integrated over the age of the solar system until the present epoch, although for this study, only a fraction of it (200 Myr) was considered, as it proved sufficient for the formation of the KB/SD and the OC reservoirs. A modified version of the *swift_rmvss4* code, part of the Swift *N*-body integration package (H. F. Levison & M. J. Duncan 1994), was used to track the orbital evolution of planets and planetesimals. In addition to the nominal output frequency of the Swift integrator (1 Myr), the orbital elements of all planetesimals were recorded with a cadence of 1000 yr, increased to 100 yr whenever a planetesimal was detected at heliocentric distances within 23 au, where increased thermal processing was anticipated.

The second set of simulations was integrated for 1 Gyr after the instability using a modified version of the Mercury hybrid integrator (J. E. Chambers 1999). The nominal output frequency for the orbital elements of planetesimals was set at 1 Myr. Additional work to increase this cadence was done using a track record of close gravitational encounters between the planetesimals and the giant planets and a track record of passages at distances within 2.5 au of the Sun. To do so, we assumed that the semimajor axis of the planetesimals' orbit remained stable between close encounters (A. Milani et al. 1987) and that any orbital changes were due to small modifications in the eccentricity. This allowed us to reconstruct individual pathways leading to ejection with a denser, yet nonconstant, output frequency that provided a better time resolution.

### 2.1. The KB/SD Sample

The KB/SD sample consists of 383 planetesimals. They were selected at the end of the simulation, which corresponds to the current epoch, by comparing the simulated orbital elements to the actual orbital distribution of observed JFCs, obtained from the JPL Small-Body Database Search Engine.<fnref n="6" /> This choice ensured that the selected planetesimals would eventually return to the inner solar system as “active” JFCs, with perihelia within 2.5 au, at least once in their lifetime and a known total absolute magnitude ($H_t$) smaller than 10.9, two conditions that allowed a better characterization of the sample.

Figure 1 shows the distribution of the KB/SD sample on the semimajor axis-eccentricity plane at the beginning ($t=0$ yr, i.e., the timing of the giant planet instability) and what we considered to be the end of the simulation in this work ($t=200$ Myr). We note that the largest fraction of planetesimals (~47%) populates the inner scattered disk, defined by orbits with $50 \text{ au} < a < 200 \text{ au}$. A fraction of 41% of the planetesimals ends up in the classical Kuiper Belt ($a < 50$ au), and a small percentage (2%) ends up at the outer scattered disk ($200 \text{ au} < a < 1000 \text{ au}$). Interestingly, a nonnegligible number of planetesimals (~9%) ends up in the OC ($a > 1000 \text{ au}$) on low-inclination orbits

<footnote n="6">⁶ http://ssd.jpl.nasa.gov/sbdb_query.cgi</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Orbital distribution in the semimajor axis and eccentricity plane of the KB/SD (black dots) and OC (red triangles) planetesimals at the beginning (t = 0 yr; top panel) and the end of the simulation (t = 200 Myr; bottom panel).</figcaption></figure>

($<$1°0), from which they are able to return into short-period and low-inclination orbits (V. V. Emel'yanenko et al. 2013; D. Nesvorný et al. 2017).

In the upper panel of Figure 2, a typical example of a planetesimal’s implantation in the inner scattered disk is presented. Departing from a low-eccentricity ($e \simeq 0.1$) orbit with a semimajor axis of $\sim$28 au, this planetesimal evolved dynamically under the influence of Neptune. Through a series of close encounters, its orbital elements were excited until the planetesimal reached the inner parts of the scattered disk

($a \simeq 90–100$ au), where it settled to high-eccentricity orbits ($e \simeq 0.5–0.7$) after the first $\sim$25 Myr of the simulation.

## 2.2. OC Sample

The OC sample consists of 360 planetesimals. It was calibrated with a sample of LPCs so that (1) the sample was of comparable size to the KB/SD sample, and (2) upon their return to the inner solar system, all its members had minimum perihelion distances within 10 au, with a linear cumulative distribution function. These arbitrary conditions were chosen
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example of the planetesimals' implantation on the Kuiper Belt (upper panel) and the OC (lower panel) on the semimajor axis vs. eccentricity plane. The color code provides the timescale of the simulation. The orange clusters denote the orbital evolution of Jupiter, Saturn, Uranus, and Neptune.</figcaption></figure>

due to the lack of a comprehensive survey for the orbital distribution of LPCs (D. Vokrouhlický et al. 2019). Figure 1 shows the initial and final distribution of the OC sample as well. Despite the relatively small number of planetesimals, we can observe the formation of the OC after 200 Myr. At that moment, the majority of particles (61%) are already placed in the OC ($a > 1000$ au), populating for the most part (75%) the region beyond the Oort peak at 10,000 au. Of the remaining 39%, some planetesimals are observed to be on transient orbits with $a > 600$ au and high eccentricities, while some populate the region with $a < 1000$ au, from where they are expected to return to the inner solar system as Halley-type comets.

The lower panel of Figure 2 shows a typical example of orbital evolution for a planetesimal of this sample as it is implanted in the OC. It started on a low-eccentricity ($e \approx 0.08$) orbit with a semimajor axis of ~28 au, and then, after a series of encounters with Neptune and Uranus, this planetesimal was

scattered inward and came under the dynamical control of Saturn after ~50 Myr. Through a series of close encounters with Saturn, it was eventually scattered toward the OC rather quickly, ~60 Myr after the start of the simulation. It reached an orbit beyond the Oort peak before it settled in the inner OC ($a < 10,000$ au) on a high-eccentricity orbit.

## 2.3. Ejected Planetesimals

The sample of ejected planetesimals consists of 675 particles. These planetesimals interacted with the giant planets after the disk was dispersed and were eventually ejected out of the solar system without ever reaching one of the reservoirs. A planetesimal was considered to be ejected when it reached distances larger than 1 pc (206, 265 au). Although this simulation was integrated for 1 Gyr, the ejection timescales were in general significantly shorter and varied between
particles, as each planetesimal's trajectory is unique and determined by the stochastic nature of its encounters with the giant planets. On average, a planetesimal was ejected after ~32 Myr (the median time was ~17 Myr) from the beginning of the simulation. The first recorded ejection took place at ~1 Myr and the last at ~186 Myr.

### 3. Dedicated Thermal and Ice Evolution Models

Our goal was to evaluate the degree of thermal processing of planetesimals during the phase of implantation or their ejection from the solar system. To this end, a dedicated thermal evolution model was designed, adapted to the needs of the coupling to the trajectories proposed by the N-body simulations. This thermal evolution model, for reasons that will become obvious in the current section, focuses mainly on two main processes of the thermally induced evolution: the energy and gas diffusion inside a porous nucleus. In this section, we describe the key elements and assumptions of our model and the numerical limitations that led us to the latter.

#### 3.1. Thermal Evolution Model

To describe the heat diffusion inside the nucleus of the simulated planetesimals, our thermal evolution model solves the energy conservation equation,

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

where $\rho_{\text{bulk}}$ (kg m⁻³) is the nucleus's bulk density, c (J kg⁻¹ K⁻¹) is the nucleus's material heat capacity, T (K) is the temperature, and $\kappa$ (W K⁻¹ m⁻¹) is the nucleus's effective thermal conductivity. The term S corresponds to the collective contributions of internal heat sources and sinks related to a variety of physical processes such as the latent heat consumed or released during the sublimation and condensation of ices (e.g., J. Benkhoff & T. Spohn 1991), heat transport by gas diffusion in the porous matrix (e.g., J. Benkhoff & W. F. Huebner 1995), internal heating by radioactive elements (e.g., D. Prialnik et al. 1987), or energy released during the exothermic process of crystallization of amorphous water ice (e.g., B. Schmitt et al. 1989).

The planetesimals were modeled as spherical, porous aggregates of dust and ices, with a radius of 5 km. The porosity was fixed to a constant value of 75%, based on porosity estimates for 67P/Churyumov-Gerasimenko (H. Sierks et al. 2015; L. Jorda et al. 2016). Thermophysical parameters were chosen as typical average values from the published literature (e.g., W. F. Huebner et al. 2006). For the dust component, a mixture of refractory and organic materials was considered. For the icy component, a selection of the more abundant ice species detected on comets (e.g., D. Bockelée-Morvan & N. Biver 2017) was selected: CO, which was also considered representative of hypervolatile species; CO₂, representative of moderate volatile species; and H₂O ice in—at least initially—amorphous form. Given the broad range of values and uncertainties regarding the abundances of those species on cometary populations (e.g., T. Ootsubo et al. 2012; O. Harrington Pinto et al. 2022) or individual comets (e.g., M. Choukroun et al. 2020; N. Biver et al. 2021), a unique, simplified initial composition was chosen. The refractories-to-ice mass ratio was set to 1, a slightly lower value than the predictions for ice-rich planetesimals at the

beginning of the simulations (e.g., M. Choukroun et al. 2020). The relative abundances to H₂O in the model were preset as follows: CO/H₂O = 5% and CO₂/H₂O = 15%, implying an initial composition of 42% H₂O ice, 6% CO₂ ice, and 2% CO ice.

The terms appearing on the left-hand side of Equation (1) ($\rho_{\text{bulk}}$, c, and the conductivity $\kappa$ of the solid component) were modeled as mass-weighted averages of the individual values of their respective components as those were presented in the previous paragraph. A selection from the most up-to-date values available in the published literature was used wherever possible. This selection is summarized in Table 2 of the Appendix. Appropriate corrections to account for the porous nature of the nucleus were introduced following typical considerations for the modeling of cometary nuclei (D. Prialnik et al. 2004; W. F. Huebner et al. 2006).

Of these terms, the most crucial in our thermal evolution model was the effective thermal conductivity, which controlled the amount of energy transferred toward a planetesimal's interior. A variety of theoretical approximations for the calculation of the effective thermal conductivity exist in the published literature. For a comprehensive review of the different approximations and their historical evolution, we refer the reader to D. Prialnik et al. (2004) and W. F. Huebner et al. (2006). In this work, we used Russel's formula (H. W. Russell 1935) as a first reduction factor $\phi$ to account for the porous nature of cometary material (see also the Appendix). Considering that the solid matrix is composed of small icy dust grains, an additional reduction known as the Hertz factor $h$ was introduced to describe conduction only through a small fraction of the grains' total surface area. This correction is calculated as the ratio of the contact area to the cross section of a spherical grain (K. J. Kossacki et al. 1999). Laboratory measurements give a wide variety of values, ranging from 10⁻¹ to 10⁻⁴ (W. F. Huebner et al. 2006). However, a consensus seems to have been established around values of the order of 10⁻² (e.g., W. F. Huebner et al. 2006; B. Gundlach & J. Blum 2012), which is the value used in the model.

Accounting for the above corrections and the contribution of the radiative transport through the pores given in Equation (A2), the final expression for the effective thermal conductivity is

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

The conductive term in Equation (2) dominates over the radiative term at low temperatures ($T < 100$ K). Above this threshold, the contribution of the radiative term becomes more significant and dominates above 150 K, as the heat is mainly transported by radiation through the pores (A. Gkotsinas et al. 2022). For a porosity of 75% and a Hertz factor $h = 10^{-2}$, our resulting thermal conductivity falls in the 10⁻³–10⁻⁴ W K⁻¹ m⁻¹ range, in good agreement with estimates for comet 67P/Churyumov-Gerasimenko (H. U. Keller et al. 2015).

#### 3.2. Coupling Thermal Calculations to N-body Simulations

The coupling between a planetesimal's orbital and thermal evolution is made through the energy balance at the surface, i.e., the surface boundary condition of the thermal
evolution model,

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

where the insolation on the left-hand side of the equation is a function of the Bond albedo (A), the solar constant (L<sub>☉</sub>), the heliocentric distance (d<sub>H</sub> in au), and the local zenith angle (ζ). We note that the solar constant is reduced to a percentage of its current value (0.8 L<sub>☉</sub>), where L<sub>☉</sub> = 1360 W m<sup>-2</sup>, in a coarse attempt to account in a simplified manner for the changes of the solar luminosity in the early stages of our solar system (D. O. Gough 1981; N. Lagarde et al. 2012). More elaborate approximations exist and have been employed on other models (see, for example, B. J. R. Davidsson 2021), but we consider that for a first-order assessment of the early thermal evolution, this simplified value is sufficient. Moreover, the local zenith angle is set to 0, meaning that our simulations focus on the processing around the subsolar point. This choice implies that we are focusing on the maximum effect of the solar illumination, as it is quantified for the purpose of this work (see Section 3.2). This end-member scenario is not necessarily characteristic of other areas of cometary nuclei, such as the poles, where lower temperatures, and as a consequence less processing, are expected. The thermal emission, on the right-hand side of the equation, is a function of the emissivity (ε), the Stefan-Bolzmann constant (σ), the temperature (T in K), and the heat flux toward the interior of the thermal conductivity (κ in W m<sup>-1</sup> K<sup>-1</sup>) at the surface.

More specifically, the coupling is realized through an accurate calculation of the heliocentric distance (d<sub>H</sub>) at each and every time step of the dynamical simulations, appearing in the calculation of the received solar energy. Obviously, such a calculation is not straightforward, as information on the position of planetesimals is only available at specific intervals.

In long-term simulations, it is a common practice to replace eccentric orbits with circular ones, assuring the conservation of a key quantity, such as the total amount of solar energy received over one orbital period (e.g., D. Prialnik & E. D. Rosenberg 2009; A. Guilbert-Lepoutre 2012; A. Gkotsinas et al. 2022; E. Lilly et al. 2024). Although such an approximation results in a loss of information, mainly related to perihelion passages and seasonal activity, it provides not only a nonnegligible improvement in numerical efficiency but also a very elegant solution to a variety of issues raised by the use of eccentric orbits. For example, it lifts the ambiguity introduced by the lack of information on the direction of the comet nucleus, i.e., inbound versus outbound, during a given calculation time step (T. K. Safrit et al. 2021). In this work, following a previous study (A. Gkotsinas et al. 2023), we employ the effective thermal distance r<sub>T</sub>, an expression proposed by A. Méndez & E. G. Rivera-Valentín (2017), following previous similar calculations for the annual mean insolation and the effective surface temperature (e.g., W. R. Ward 1974), producing circular orbits with the same average equilibrium temperature for an eccentric orbit with a given semimajor axis (a) and eccentricity (e),

$$ r_T = a \left[ \frac{2\sqrt{1+e}}{\pi} E \sqrt{\frac{2e}{1+e}} \right]^{-2} \quad (4) $$

$$ \approx a \left( 1 + \frac{1}{8} e^2 + \frac{21}{512} e^4 + \mathcal{O}(e^6) \right), \quad (5) $$

with E, the complete elliptic integral of the second order. As explained, this approximation does not allow the resolution of seasonal variations of the insolation related to eccentricity or the consideration of seasonal variations related to obliquity. The same goes for diurnal variations, which are unique to each planetesimal and cannot be studied in a statistically significant manner.

### 3.3. Internal Energy Sources and Sinks

While the orbital elements in the N-body simulations are recorded with a relatively large cadence (e.g., 100–1000 yr in the first set of simulations), thermal processes usually operate on much shorter timescales.

In Figure 3, these timescales are plotted as a function of temperature; for a hypothetical subsurface layer of thickness Δr = 1 m; with a porosity of 50%; with typical cometary values for the bulk density, the heat capacity, and the effective thermal conductivity; and for different compositions. We can easily observe the wide range of timescales on which each process operates, accelerating significantly as the temperature rises, and reducing the timescale by several orders of magnitude. In addition, we observe important differences when these processes are compared to each other at a given temperature. A striking example is the comparison between the thermal and the gas diffusion timescales, the former varying between a few months (10<sup>6</sup> s) and a little over a year, depending on the considered composition, while the latter takes place in a matter of hours or minutes, depending on the size of the pores. Sublimation is an exponential function of the temperature: CO<sub>2</sub> and H<sub>2</sub>O ices sublimate in a few minutes at ~80 K and ~180 K, respectively, but are practically insignificant below these threshold temperatures. The time needed for the crystallization of amorphous water ice at temperatures below ~80 K exceeds the solar system's lifetime but can take place in just a few years for temperatures between 100 and 110 K and in only a few hours for temperatures between 125 and 150 K (D. Jewitt 2009; D. Prialnik & D. Jewitt 2022). Accounting for processes with such diverse timescales has already been one of the main challenges of thermal evolution models (e.g., R. Orosei et al. 1999; D. Prialnik & E. D. Rosenberg 2009), even when the integration period was short and only a single object was studied. In this work, as we dealt with long-term simulations requiring a substantial time of integration and as we used a large number of objects to get a statistically significant image of the processing, further simplifying assumptions were required.

We fixed the time step for thermal calculations at 1 yr, unless smaller values were available in the case of the ejected planetesimals. This means that the thermal evolution model created smaller time steps within each N-body simulation's output that allowed a more accurate calculation of the internal temperature of the planetesimals. This value was considered a good compromise between the different timescales involved in this study. It allowed the resolution of the heat diffusion in the interior of the planetesimals, a process operating on timescales of a little over 1 yr for a planetesimal with properties close to those of comet 67P (thick red line in Figure 3) while being comparable to the N-body simulations' output cadence. On the downside, the other thermal processes could not be resolved, and they were not
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Characteristic timescales of the main physical processes taking place in a comet’s nucleus: crystallization of amorphous H₂O ice, denoted by the green dashed-dotted line; sublimation of different volatile species, denoted by the blue dashed lines; gas diffusion within the pores for pore radii between 10⁻³ and 10⁻⁶, denoted by the orange dotted lines and the orange shaded area; and heat diffusion in a nucleus of different compositions: pure dust, pure crystalline H₂O ice, and a nucleus with characteristics resembling those of comet 67P/Churyumov–Gerasimenko. The black loosely dotted lines serve as timescale indicators. Figure adapted from D. Prialnik et al. (2004) and W. F. Huebner et al. (2006).</figcaption></figure>

taken into account in our thermal evolution model. The possibility of heating due to the decay of radioactive elements was not considered as well, as on the one hand, the presence of radioactive elements remains uncertain (e.g., A.-C. Levasseur-Regourd et al. 2018; U. Malamud et al. 2022), and on the other hand, the timing of the instability, and therefore the inclusion of radioactive elements on planetesimals, is still a matter of debate (e.g., S. N. Raymond et al. 2020; G. J. Golabek & M. Jutzi 2021).

To summarize, the current configuration of our thermal evolution model neglects phase transitions, at least in a classic direct way, implying that the term S on the right-hand side of Equation (1) is 0. As a consequence, it does not faithfully describe periods of cometary activity, when thermally induced phase transitions contribute significantly to the total energy budget. Ideally, this should be obtained by the simultaneous resolution of a mass conservation equation (e.g., B. J. R. Davidsson 2021), the results of which are used in the sources and sinks term of the energy conservation equation (Equation (1)). Given the long time steps of our simulations, this classic resolution scheme was not evident, so we opted for an alternative approach to get an estimation of the ice content evolution. This approach consisted of the asynchronous coupling of an ice evolution model to the thermal evolution model, a solution initially proposed by N. Schorghofer (2008) and further developed in N. Schorghofer (2010, 2016) and N. Schorghofer & H. H. Hsieh (2018).

### 3.4. Asynchronous Ice Evolution Model

The basic idea behind the ice evolution model is that we were dealing with temperature-driven processes. Therefore, by calculating an average of the temperature variations over a given time interval (for example, a time step of 100 or 1000 yr)

can give us an estimation of the evolution of the ice content over this interval by calculating the ensuing average vapor pressures developing in a planetesimal’s interior. In the following paragraphs, a brief overview of the ice evolution model is provided. We refer the reader to the aforementioned publications for a more detailed description.

Assuming an ice interface buried at some depth at a radial distance from the center of a planetesimal’s nucleus ($r_i$), following N. Schorghofer & H. H. Hsieh (2018), we write the vapor flux at the interface as

$$J(r_i) = \rho_i \frac{dr_i}{dt}, \qquad (6)$$

where $\rho_i$ (kg m⁻³) is the ice density.

Supposing that the vapor flux is proportional to 1/r², we write the vapor flux at a random radial distance from the center as

$$J(r) = -D \frac{d\rho_v(r)}{dt} = \left(\frac{r_i}{r}\right)^2 J(r_i), \qquad (7)$$

where D (m² s⁻¹) is the vapor diffusion coefficient and $\rho_v$ (kg m⁻³) is the vapor density. The vapor diffusion coefficient can be expressed as

$$D = \frac{\pi \psi v_{th}}{8 + \pi (1 - \psi) \xi} r_p, \qquad (8)$$

with $\psi$ the nucleus’s porosity, $v_{th} = \sqrt{8k_bT/\pi m}$ (m s⁻¹) the mean thermal speed, $r_p$ (m) the pores’ radius, and $\xi$ the tortuosity, an intrinsic property of the porous medium defined as the ratio between the actual flow path to the straight distance between the two ends of this flow path (J. Bear 1988). This term seeks to
describe the “twisting” of the pores in a porous medium structure and here is set to 1, implying that the pores, or the flow path, is considered to be a straight line, a very common assumption in thermal evolution models (W. F. Huebner et al. 2006). It should be noted that other expressions for the vapor diffusion coefficient exist, such as Knudsen’s formula (S. Espinasse et al. 1991) or the Clausing equation (G. Steiner 1990), but as A. Guilbert-Lepoutre (2014) reported, the calculated vapor fluxes differ by a factor that does not exceed the order of unity.

Resolving Equation (7) for the vapor density and subsequently deriving it at $r = r_i$, supposing that the vapor pressure at the surface is 0 and at the interface equal to the saturation vapor density ($\rho_s$), we obtain the following expression for the vapor flux at the ice interface (N. Schörghofer & H. H. Hsieh 2018):

$$J(r_i) = \frac{D\rho_s}{r_i(1 - r_i/R)} \qquad (9)$$

Combining Equations (6) and (9), we can obtain an expression for the retreat rate of the ice interface in the form of a differential equation that can be integrated to a time-stepping scheme,

$$\frac{dr_i}{dt} = D \frac{\rho_s}{\rho_i r_i (1 - r_i/R)}, \qquad (10)$$

or, more conveniently,

$$r_i \left(1 - \frac{r_i}{R}\right) dr_i = D \frac{\rho_s}{\rho_i} dt. \qquad (11)$$

We note that the saturation vapor pressure ($\rho_s$) depends only on the temperature and can be calculated from the saturation vapor pressure ($P_{sat}$ (Pa)) using the average temperatures provided by the thermal evolution model at each time step provided by the dynamical simulations (100 or 1000 yr for KB/SD and OC objects and 1 Myr or smaller for ejected planetesimals).

Here, we calculate the saturation vapor pressure using a simple empirical expression obtained from the Clausius–Clapeyron equation for the CO₂ and the H₂O ice molecules (F. P. Fanale & J. R. Salvail 1984, 1990),

$$P_{sat}(T) = Ae^{-B/T}, \qquad (12)$$

where the coefficients A and B are $107.9 \times 10^{10}$ Pa and 3148 K for CO₂ and $356 \times 10^{10}$ Pa and 6141.667 K for H₂O, respectively. This expression, although relatively simple, is in excellent agreement with more complex formulae, such as the ones proposed by G. N. J. Brown & W. T. Ziegler (1980) and R. Orosei et al. (1999). However, for CO₂, this expression, using coefficients from the same work (F. P. Fanale & J. R. Salvail 1984), presents important discrepancies with more complex ones (G. N. J. Brown & W. T. Ziegler 1980; N. Fray & B. Schmitt 2009) and with the results of recent laboratory measurements (W. Grundy et al. 2024) that measured significantly lower values, by at least 2 orders of magnitude, for temperatures below 30 K. In this study, an expression proposed by G. N. J. Brown & W. T. Ziegler (1980) is employed, as it is considered a good compromise between the more complex expressions and the aforementioned laboratory

measurements,

$$\ln(P_{sat}) = A_0 + \sum_{i=1}^{5} A_i T^i, \qquad (13)$$

where the values for the coefficients $A_0$ and $A_i$ are given in Table 3 of the Appendix.

Obviously, this analytical method does not explicitly account for internal processes such as the gas flow or the recondensation below and above the ice interfaces. It does, however, remain a powerful tool for a first-order estimation of a planetesimal’s evolution, as it still takes into account the two dominant processes of the thermally induced activity: the heat and the gas diffusion inside the nucleus.

### 3.5. Supplementary Ending Conditions

As we mentioned earlier, for this work, the orbital evolution of the planetesimals was recorded for a period of 200 Myr or until the ejection of a planetesimal from the solar system. For the thermal processing simulations, we added three supplementary termination conditions to further optimize the numerical efficiency.

The first condition was introduced to terminate a simulation when no thermal modification was recorded after a number of iterations. A threshold of 100 iterations was defined, after which, if no modification of the ice interfaces was recorded, the simulation ended. This condition targeted planetesimals that quickly consumed their hypervolatile content or that were rapidly placed in faraway, relatively stable orbits, where the temperatures were too low for any further activity to take place.

The second condition forced the termination of the thermal simulations if a planetesimal’s semimajor axis exceeded 100 au. This condition concerned planetesimals that managed to “maintain” some minor activity, mainly due to small changes in the float precision, “escaping,” in a way, the first condition. As these changes concern mostly planetesimals orbiting at large heliocentric distances (>70 au), the 100 au limit was considered a safe choice. An estimation of a blackbody’s temperature at 100 au for a luminosity value of 0.8 $L_\odot$ yields ≈22 K, well below the sublimation temperatures for CO₂ and H₂O and at the limit for CO₂, which, however, was never present close to the surface at these distances. We verified that no further orbital changes significantly reducing the semimajor axis of the planetesimals occurred after the termination of our thermal simulations. Extending this limit to 200 au did not provoke any substantial change in the outcome of our simulations.

Finally, a termination condition targeting the ejected planetesimals was added, ending the thermal evolution simulations when their orbits became hyperbolic, flagged by eccentricity values above 1. This choice was made first and foremost for consistency with our orbital averaging technique, valid only for bound, eccentric orbits. Either way, from a thermal evolution point of view, it was a “harmless” condition, since the transition from eccentric to hyperbolic orbits generally occurs at large heliocentric distances where no significant thermal alteration was recorded.

The application of these ending conditions, as expected, reduced the time of the thermal evolution simulations: the evolution of KB/SD planetesimals was simulated for ~63 Myr on average (median time of ~20 Myr), with the shortest simulation lasting ~1.4 Myr. For the OC sample, the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Thermal processing of a planetesimal scattered toward the OC. Upper panel: evolution of the effective thermal radius ($r_T$). Lower panel: temperature distribution (denoted by the color code) in the planetesimal's interior at a subsurface layer of ~1000 m. The dashed, dotted, and dashed-dotted white lines show the CO, $CO_2$, and crystallization interfaces, respectively.</figcaption></figure>

corresponding numbers were considerably lower, ~23 Myr and ~11 Myr for the mean and the median simulation timescales, respectively, and ~640 kyr for the shortest simulation. These numbers provide a first hint of the differences between the two populations: planetesimals scattered toward the OC were not only scattered further out but also much faster. The effect was similar on the ejected planetesimals: a significant decrease on the thermal simulation timescale. On average, a planetesimal reached a hyperbolic orbit in ~15 Myr, with the median value reduced to ~10 Myr. The shortest ejection was recorded to take place in just ~1000 yr.

# 4. Results

## 4.1. Individual Planetesimal Example

Figure 4 shows an example of the coupling between the thermal and ice evolution models and the orbital evolution of a planetesimal scattered toward the OC. We plot data relative to a 1.4 km subsurface layer only, allowing a better visualization of the temperature and ice interfaces' evolution. The thermal evolution of this planetesimal was computed for a period of ~1 Myr, before the planetesimal was scattered beyond the threshold of 100 au. The lower panel of Figure 4 displays the

temperature distribution in the planetesimal's interior, calculated by the thermal evolution model. Plotted over this temperature distribution are the evolution of the positions of the three ice interfaces, as calculated by the ice evolution model.

After its formation, the planetesimal from Figure 4 was scattered inward by the giant planets, being thrust interior to Jupiter's orbit at $t \approx 600$ kyr before being scattered outward starting after $t \approx 800$ kyr and eventually being implanted in the OC. The time that it spent under the dynamical control of Jupiter was an intense heating period that lasted ~200 kyr, characterized by a very low effective thermal radius. Five short-period approaches can be easily identified during this period (upper panel of Figure 4), associated with intense heating episodes in the lower panel of Figure 4. While CO appears to sublimate continuously within 30 au, these approaches provoked the sublimation of $CO_2$ ice and the crystallization of amorphous water ice from a substantial subsurface layer of ~90 m and ~210 m, respectively. It is important to clarify that these drops in the effective thermal radius and the ensuing heating episodes should not to be confused with the increased heating sustained during perihelion passages. The drops in effective thermal radius represent the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Cumulative distribution of the CO ice interface depths for the members of the three studied populations, KB/SD (red dashed line), OC (green dashed–dotted line), and ejected (black dotted line) planetesimals, at the end of their reservoir implantation and ejection phases, respectively.</figcaption></figure>

planetesimal's tendency to orbit closer to the Sun over periods of time that are much longer than a typical cometary orbital period, comparable to the *N*-body simulation time steps of 100 and 1000 yr. As the scattering toward the OC began, the planetesimal started to cool down, putting an end to the progression of the CO₂ and amorphous H₂O ice interfaces. Following the heat diffusion in deeper parts of its interior, the CO interface continued to regress for some time before becoming almost stagnant as the planetesimal moved beyond 100 au, and the simulation ended.

## 4.2. Processing of Populations

Having presented an example of the coupling between the evolution models and the *N*-body simulations, we now focus our attention on the three studied populations, searching for statistically significant differences that may characterize their degree of processing. The level of processing of each planetesimal, and therefore of each population, was tracked by the evolution of the ice content, as described by the retreat of the ice interfaces.

### 4.2.1. Hypervolatile Content

Figure 5 shows the cumulative distribution of the depths of the CO ice interface at the end of the thermal evolution simulations for all the planetesimals of the ejected, OC, and KB/SD populations. A clear distinction in the degree of processing for the three populations is observed. Even though the ejected population had significantly shorter lifetimes in the solar system, it is by far the most processed population. In ~76% of the ejected planetesimals, the CO ice interface had already retreated to the center of their nuclei before their ejection from the solar system. The survival of pure CO ice close to the surface proved to be a highly unlikely event, with

only ~11% of the planetesimals preserving some CO within 2500 m from the surface.

Planetesimals implanted in the KB/SD were also substantially processed. In ~51% of the planetesimals, the CO interface receded to the center of their nuclei at the end of the simulations. However, ~33% of the KB/SD planetesimals retained some quantity of their initial pure CO ice in the upper half of their nucleus.

The least processed population, and the one with the most chances of preserving an important amount of its pure condensed hypervolatile content, was the OC population, the progenitors of LPCs. Although 38% of the planetesimals were unable to retain their initial CO ice content, the chances of its survival were significantly higher than for the other populations. The possibility of finding the surviving pure CO ice close to the surface was higher, as ~46% of the OC planetesimals maintained the CO ice interface in the upper part of their nuclei (within 2500 m from the surface) at the end of the simulations. We should clarify at this point that this finding does not imply that any future activity would be triggered from such depths. It simply shows that an important reserve of pure CO ice can survive closer to the surface when compared to the other populations, and it can become an important activity driver later on, as it can gradually reach the surface, provided the necessary temperature and pressure conditions. In addition, one should take into account the particularities of the return of these objects to the inner solar system, which is quite different between the two populations —slow in the case of KB/SD objects versus sharp in the case of OC objects—and is expected to increase these differences between the two populations (e.g., A. Gkotsinas et al. 2022; N. A. Kaib 2022).

> In addition, we should recall at this point that our thermal evolution model did not explicitly account for a number of internal processes by simultaneously resolving the heat and gas
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Cumulative distribution of the CO₂ ice interface depths for the members of the three studied populations, KB/SD (red dashed line), OC (green dashed-dotted line), and ejected (black dotted line) planetesimals, at the end of their reservoir implantation and ejection phases, respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Cumulative distribution of the amorphous water-ice crystallization interface depths for the members of the three studied populations, KB/SD (red dashed line), OC (green dashed-dotted line), and ejected (black dotted line) planetesimals, at the end of their reservoir implantation and ejection phases, respectively.</figcaption></figure>

diffusion equations (as done by B. J. R. Davidsson 2021 and A. Parhi & D. Prialnik 2023, for example). Instead, it relied on a simplified thermal evolution calculation coupled asynchronously to an ice evolution model. The latter assumed that sublimated gases were vented at the surface, similarly to some previous studies (Y.-J. Choi et al. 2002; J. K. Steckloff et al. 2021). However, recondensation both above and below an ice interface is expected if the appropriate temperature and

pressure conditions are met (see A. Parhi & D. Prialnik 2023, in particular). With that in mind, we must clarify that the term “loss” for a given volatile does not necessarily signify its complete depletion. It would be more precise to say that it describes a modification of the initial mass fraction because of thermal processing. Additional seasonal or diurnal cycles, acting on shorter timescales, should allow for the local depletion or enrichment of subsurface layers.
<table><caption>**Table 1**
Final Depths of the Examined Ice Interfaces—CO, CO₂, and Amorphous H₂O Ice—at the End of the Simulations for the Three Populations: KB/SD, OC, and Ejected Planetesimals</caption><thead><tr><th>Ice Species</th><th>Max.<br>(m)</th><th>Min.<br>(m)</th><th>Med.<br>(m)</th><th>Mean<br>(m)</th><th>σ<br>(m)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">KB/SD</td></tr><tr><td>CO</td><td>5000.0</td><td>392.5</td><td>5000.0</td><td>3577.9</td><td>1619.6</td></tr><tr><td>CO<sub>2</sub></td><td>65.4</td><td>0.5</td><td>0.5</td><td>0.7</td><td>3.3</td></tr><tr><td>H<sub>2</sub>O<sub>am</sub></td><td>184.5</td><td>0.3</td><td>0.3</td><td>0.8</td><td>9.5</td></tr><tr><td colspan="6" style="text-align: center;">OC</td></tr><tr><td>CO</td><td>5000.0</td><td>287.1</td><td>2816.7</td><td>3081.9</td><td>1688.3</td></tr><tr><td>CO<sub>2</sub></td><td>88.5</td><td>0.5</td><td>0.5</td><td>2.5</td><td>9.7</td></tr><tr><td>H<sub>2</sub>O<sub>am</sub></td><td>232.5</td><td>0.3</td><td>0.3</td><td>6.1</td><td>27.3</td></tr><tr><td colspan="6" style="text-align: center;">Ejected</td></tr><tr><td>CO</td><td>5000.0</td><td>167.4</td><td>5000.0</td><td>4420.6</td><td>1219.2</td></tr><tr><td>CO<sub>2</sub></td><td>5000.0</td><td>0.5</td><td>0.5</td><td>45.2</td><td>430.9</td></tr><tr><td>H<sub>2</sub>O<sub>am</sub></td><td>1052.5</td><td>0.3</td><td>0.3</td><td>16.5</td><td>66.6</td></tr></tbody></table>

### 4.2.2. Moderate- and Low-volatility Content

Figures 6 and 7 show the cumulative distribution of the depths of CO₂ and amorphous H₂O ice at the end of the thermal evolution simulations, similar to Figure 5. Compared to alterations on the hypervolatile content, here, only a small fraction of planetesimals from each population undergoes the high-temperature processing leading to the sublimation of CO₂ ice and the crystallization of amorphous H₂O ice. Nonetheless, in terms of CO₂ and amorphous H₂O ice, ejected planetesimals are again more processed on average than the “surviving” planetesimals implanted in the KB/SD and the OC. A total of 18% of the ejected planetesimals lost part of their CO₂ ice content, and in ~25%, the amorphous H₂O ice crystallized in a subsurface layer that did not exceed ~100 m. Only a tiny fraction (~2%) of ejected planetesimals underwent extensive processing, with ice interfaces receding as deep as ~1000 m.

The affected planetesimals in the KB/SD and the OC populations were significantly lower: only ~11% and ~9% in the case of the OC (for CO₂ and amorphous H₂O ice, respectively) and just ~1.5% for the KB/SD planetesimals. These reduced fractions of affected planetesimals in the OC and KB/SD populations are a natural consequence, on the one hand, of the higher temperatures needed for these processes to take place and, on the other hand, of the scattering procedure that led them to the source reservoirs, as we briefly showed in Figure 2 and as we explain in the next paragraphs. An overview of our simulations' results can be found in Table 1.

### 4.2.3. The Imprint of the Orbital Evolution on the Thermal Processing

As shown in Figure 2, KB/SD planetesimals were typically scattered outward by Neptune and Uranus (e.g., D. Nesvorný 2015). This happens over relatively long timescales, spending considerable time in the Kuiper Belt and the inner parts of the scattered disk (with <100 au), where temperatures are sufficient to trigger the sublimation of CO. However, this same pattern of orbital evolution resulted in the depletion of their condensed hypervolatile content that protected their moderate- and low-volatility species. As they were mostly under the influence of the migrating Neptune, these

planetesimals rarely roamed the area within Neptune's or Uranus's orbits, where temperatures are high enough for the sublimation of CO₂ ice and the crystallization of amorphous H₂O ice, before their implantation on stable orbits in the KB/SD.

Planetesimals implanted in the OC were mainly scattered outward by Jupiter and Saturn (e.g., L. Dones et al. 2015; N. A. Kaib & K. Volk 2022). This difference implies more passages in the inner parts of the giant planet area and orbits closer to the Sun than the KB/SD population. These passages can eventually lead to intense, generally short-lived, heating episodes that can affect species of lower volatility, affecting a shallow subsurface layer, as in the example of Figure 4. This also explains why, if the CO₂ and crystallization interfaces were used as processing indicators, the KB/SD population appears to be less processed than the OC one, unlike the estimation based on the processing of CO.

The same pattern seems to apply for the ejected planetesimals, which were for the most part scattered outward by Jupiter. Although their lifetimes in the solar system were much shorter than that of OC planetesimals (see Section 3.5), they spent them almost entirely in the inner parts of the solar system, within the orbit of Jupiter, where they were exposed to much higher temperatures, leading to extensive processing, affecting all the primordial pure ice inventories.

## 5. Discussion

### 5.1. On the Survival of CO Ice in the KB/SD

Evaluating the long-term global loss of condensed hypervolatile species in the progenitors of JFCs has been the subject of several studies (M. C. De Sanctis et al. 2001; Y.-J. Choi et al. 2002; B. J. R. Davidsson 2021; J. K. Steckloff et al. 2021; A. Parhi & D. Prialnik 2023). It requires models simultaneously solving both the heat and gas diffusion equations (see also Section 3) and simulations integrated on long timescales to evaluate whether this sublimation process is completed. For example, M. C. De Sanctis et al. (2001) performed such simulations on large objects (R = 80 km) evolving on stable orbits, so most likely postreservoir implantation, at distances between 39 and 43 au to assess the influence of insolation and radiogenic heating. After 10 Myr of simulation, they found that ~30% of the initial CO content was lost. Based on the exponential decline of the sublimation rate, they assumed that CO would be able to survive at the innermost parts of the objects.

Similar studies were conducted by Y.-J. Choi et al. (2002) and J. K. Steckloff et al. (2021). They found that condensed CO should be completely lost, neglecting, however, gas diffusion and assuming simply that the sublimating gas would vent at the surface. B. J. R. Davidsson (2021) simulated the coupled heat and gas diffusion on timescales long enough to actually assess the loss of pure condensed hypervolatile species. They found that JFC-sized objects (R = 4 km) should lose all CO within 0.1 Myr, unless CO is trapped in some other medium such as amorphous H₂O ice or CO₂ ice. That timescale significantly increased with the objects' size, as 70–200 km objects were found to lose all condensed CO within 10–13 Myr. A. Parhi & D. Prialnik (2023) considered in their study more complex mixtures. They showed that all JFC-sized objects should lose all condensed hypervolatile species, and that larger objects should not retain a primitive volatile mass
ratio in their subsurface layers. Their work demonstrated an interesting sublimation and recondensation pattern, where volatile depletion proceeded from internal bulk sublimation and from a sublimation front receding from the surface. A layer significantly enriched in hypervolatiles developed just below the sublimation front, which was located close the surface.

What this work adds to this picture is the orbital evolution of planetesimals. All studies mentioned above assumed icy objects on fixed orbits for the duration of their thermal simulations. Despite the different assumptions between the thermal evolution models and the considerations of the aforementioned analytical studies, our results are generally in agreement with the prediction of pure condensed CO ice depletion at heliocentric distances up to 45 au. However, beyond this distance, the chances for condensed CO ice survival increase. We suggest an increasing survival rate for condensed CO for scattered disk objects and especially those placed at distances beyond 200 au, a region considered to be the main source for objects evolving into Jupiter-crossing orbits (e.g., D. Nesvorný 2015). It is the consideration of their orbital evolution, by means of N-body simulations, that revealed the effect of the chaotic nature of the trajectories on the thermal evolution of these planetesimals, bringing them back and forth from warmer to colder regions in a random way (A. Gkotsinas et al. 2022). This proved to be a key factor in the retention of condensed CO ice in their interiors. Unlike previous studies, our objects underwent multiple heating and cooling episodes during their implantation on the outer solar system reservoirs, which acted as a significant “protection” mechanism for their buried, condensed hypervolatile content.

## 5.2. On the Processing of OC Objects

The few studies that examined the thermal processing of objects in the OC focused mainly on the effects of radioactive heating (J. Haruyama et al. 1993; D. Prialnik & M. Podolak 1995). Given the completely different set of initial assumptions, a comparison between their results and ours does not seem pertinent. A more prominent comparison can be made with some recent predictions by B. J. R. Davidsson (2021), who noted that the increased timescales for CO loss in large objects (R = 70–200 Km) imply a dynamical evolution effect on the survival of CO, especially if they are scattered quickly toward the OC. This is what our results confirm, as the OC planetesimals proved to be the least processed population, a result of a relatively quick outward scattering. Indeed, the implantation on the OC took place on different timescales and appeared to be a slow as much as a quick procedure (see Section 3.5), leading to a variety of thermal processing outcomes, among which ~60% of the OC planetesimals managed to retain a part of their primordial pure condensed CO ice. This is in disagreement with predictions for a lack of hypervolatile species in OC comets, suggested by a recent analytical study (C. M. Lisse et al. 2022). However, the differences in our results can be tracked back to dynamical timescales, as C. M. Lisse et al. (2022) assumed that the implantation of planetesimals toward the OC starts after objects spent hundreds of Myr in the KB/SD region. Given that there is strong evidence that the giant planet instability took place no later than 100 million yr after the start of solar system formation (A. Morbidelli et al. 2018; D. Nesvorný et al. 2018; S. J. Mojzsis et al. 2019), and quite possibly within less than 10 million yr (e.g., M. S. Clement et al. 2018;

R. R. de Sousa et al. 2020; B. Liu et al. 2022), this working assumption by C. M. Lisse et al. (2022) seems rather unlikely.

Our prediction for condensed CO ice survival in the OC is consistent with the growing number of observations of suspected CO-driven activity on inbound LPCs, such as C/2017 K2 (PANSTARRS) (K. J. Meech et al. 2017), C/2010 U3 (Boattini) (M.-T. Hui et al. 2019), and C/2014 UN271 (Bernardinelli-Bernstein (T. L. Farnham et al. 2021). At the same time, the possibility of CO depletion or low CO abundance ratios, as observed, for example, on comet C/2021 A1 (Leonard) (S. Faggi et al. 2023), is not excluded at all.

In addition, our results provide an explanation of the observed differences in the CO/CO₂ and CO/H₂O abundance ratios of LPCs (M. F. A'Hearn et al. 2012; O. Harrington Pinto et al. 2022), even if the only available data were obtained at short heliocentric distances (up to ~6 au) and even if those are based on production rates that cannot be directly calculated by our model. According to our simulations, these differences are an inherent consequence of the planetesimals' orbital history, an explanation that seems to be more plausible than a scenario for different formation areas for KB/SD and OC planetesimals, as suggested by M. F. A'Hearn et al. (2012).

This explanation can be extended to the production rate differences between dynamically new (DN) and dynamically old (DO) LPCs, observed initially by M. F. A'Hearn et al. (1995) and confirmed recently by O. Harrington Pinto et al. (2022). Both studies observed that DN LPCs produce more CO₂ than CO, while DO LPCs appear to be more CO-dominant. This unexpected observation can be explained by the different trajectories leading to implantation on the OC. Assuming a typical value of 10,000 au for the OC spike separating DN from DO LPCs, we can see in Figure 8 that planetesimals that are statistically more likely to visit the planetary region for the first time (DN) are more CO-processed on average than planetesimals implanted below the Oort spike, an area that is considered to be the source of DO LPCs. This result does not exclude the possibility of further processing in the OC by galactic cosmic rays, as has been previously proposed to explain this peculiar observation (O. Harrington Pinto et al. 2022), but provides an additional and more concrete explanation based on the implantation specifications on the OC.

We should note at this point that this explanation does not automatically apply to comparisons between JFCs and LPCs, as one might expect. If, for LPCs, our assumption seems more reliable given the abrupt nature of their return to the inner solar system, we cannot say the same thing for JFCs, which follow much more complicated and long trajectories toward the planetary region that is expected to provoke additional thermal processing (e.g., A. Gkotsinas et al. 2022), further affecting the chemical and physical properties of the comets. To make such comparisons, further simulations accounting for processing during the return to the inner solar system are necessary, accompanied by stronger observational constraints on CO and CO₂ production rates, which for the time being are limited at distances within 3.5 au (O. Harrington Pinto et al. 2022).

## 5.3. On Drivers for Cometary Activity

In agreement with these previous studies (see Section 5.1), no significant alterations on the moderate- (e.g., CO₂ or amorphous H₂O ices) and low-volatility (e.g., H₂O ice) content were recorded on planetesimals ending up in the KB/SD, when the only energy source considered is the solar irradiation.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Cumulative distribution of the CO ice interface depths on planetesimals implanted at distances below the Oort spike ($a < 10,000$ au) and considered to be DO (orange dotted line) and planetesimals implanted beyond the Oort spike, considered to be DN (blue dashed line).</figcaption></figure>

Alterations on the entire icy inventory, if any, should be the result of heating due to the decay of radioactive elements, affecting large objects more significantly ($R \ge 40$ km) with low effective thermal conductivity (e.g., M. C. De Sanctis et al. 2001; Y.-J. Choi et al. 2002; U. Malamud et al. 2022). We found that for most KB/SD (~98%) and OC objects (~90%), H₂O ice remained in its initial amorphous phase during their implantation on their respective reservoirs, meaning that until their return to the inner solar system as JFCs and LPCs, it is highly unlikely to undergo any further alterations. This result implies different activity drivers upon the return of these objects in the inner solar system as JFCs or LPCs. As CO seems to be more abundant and closer to the surface on OC objects, we would expect that activity, at least at large heliocentric distances, to be driven by CO sublimation. This prediction is in line with a growing number of suspected CO-driven activity on inbound LPCs (e.g., K. J. Meech et al. 2017; M.-T. Hui et al. 2019; T. L. Farnham et al. 2021). On the other hand, CO ice, although present, seems to be less abundant on KB/SD objects, while the CO₂ and amorphous H₂O ice content appeared to remain almost intact and much closer to the surface. This suggests a larger diversity of potential drivers upon the return of KB/SD objects on the inner solar system, with less pronounced CO-driven activity.

We would like to clarify here that the results presented above and the ensuing interpretations should not be viewed as an attempt to explain a vast and very diverse body of observations regarding the activity of comets. This would require a more complicated evolution model, and even then, given the large number of free or unknown parameters, different outcomes should be expected. Instead, we aim to highlight the importance of the orbital evolution, a key aspect, so far disregarded, that is intimately linked to the thermally induced physical and chemical evolution of cometary nuclei. By doing so, we showed that a variety of contradictory results, such as the predictions for CO depletion (e.g., B. J. R. Davidsson 2021;

C. M. Lisse et al. 2022; A. Parhi & D. Prialnik 2023) or the observation of long-distance CO-driven activity (e.g., K. J. Meech et al. 2017), can be explained in a very natural way. In addition, and despite the simplifications and the working assumptions, we showed that we can clearly identify statistically distinct behaviors between the different populations, even though they were formed and evolved from the same region of the protoplanetary disk.

### 5.4. A Critical Review of Our Method

In Section 3, a novel approach was proposed to examine the thermal processing of planetesimals by taking into account their long-term orbital evolution. In this approach, a number of simplifying assumptions were employed, rendering the coupling of the orbital to the thermal evolution feasible. In the next paragraphs, we would like to examine the effect of these simplifications on the outcomes of our simulations to better evaluate them.

#### 5.4.1. Implications of Orbital Averaging

One of the most prominent assumptions of our method is the use of equivalent average equilibrium temperature circular orbits. Although it proved to be a powerful tool from a numerical point of view, solving, for instance, several issues related to the output cadence of *N*-body simulations (see also A. Gkotsinas et al. 2023), their use came with an important setback: they poorly accounted for high-temperature processes such as the sublimation of H₂O ice.

As explained in A. Gkotsinas et al. (2023), all averaging techniques systematically place objects closer to their aphelion distance, resulting in the underestimation of the surface temperature, especially close to perihelion. As a consequence, the simulated surface temperatures were insufficient for the onset of H₂O ice sublimation and, in some cases, of CO₂ sublimation or the crystallization of amorphous H₂O ice.
The H₂O ice sublimation from the surface and the associated erosion should bring the buried ice interfaces closer to the surface, prompting an enhancement of the different processes as the incoming solar energy reaches these interfaces more easily. Unfortunately, such an effect is not accounted for in this study. However, given the low sublimation temperature of condensed CO, our results for this ice are not substantially affected, as the CO₂ and amorphous H₂O ice interfaces might be. With that in mind, our results for OC planetesimals that had prolonged residence times in the inner planetary region should be considered rather conservative on moderate- and low-volatility species.

### 5.4.2. Implications of the Thermal Modeling Assumptions

In Sections 3.1 and 4.2.1, we described some of the impacts that our working assumptions for the thermal evolution model may have on our results. We highlighted in particular the short-term processes (of the order of minutes or hours) that are lacking in our methods. Indeed, ice recondensation above and below an ice interface should trigger a substantial redistribution of the ice species in a cometary interior. Moreover, we pointed out that our model examines the theoretically maximum processing of an area at the subsolar point, receiving an average amount of solar energy that is decided from Equation (5), meaning that no amplitude effects from any temporal or geographical variations were considered.

In Section 3.4, we discussed in length the effect of an accurate estimation of the saturation vapor pressure on CO ice conservation. Using classical formulae from F. P. Fanale & J. R. Salvail (1990; see also D. Prialnik et al. 2004), the chances of CO survival, even in the deepest part of the nucleus, were close to zero. With more accurate estimations, such as the ones from G. N. J. Brown & W. T. Ziegler (1980), N. Fray & B. Schmitt (2009), and W. Grundy et al. (2024), we found that CO ice survival is instead possible. This observation points to the necessity of obtaining better quality constraints on parameters involved in cometary thermal evolution models (either from laboratory experiments or through observations), as slight changes in their values can radically change the outcomes of a given simulation. For example, some studies model the dust component as large aggregates of millimeter size (e.g., B. Gundlach & J. Blum 2012; U. Malamud & D. Prialnik 2015), implying also larger pores, which are usually modeled to have a radius of the same order of magnitude as the grains (e.g., W. F. Huebner et al. 2006). Here, the dust component was modeled as a powder with grain sizes of the order of 1 µm, a typical working assumption on many thermal evolution models (e.g., A. Enzian et al. 1997; W. F. Huebner et al. 2006; S. Loveless et al. 2022). We tested the impact of such an assumption on a smaller sample of KB/SD planetesimals, using larger pores of the order of millimeter size to test the other end-member scenario. As expected, the effects on the diffusivity and subsequently on the ice loss (see Equations (8) and (11)) were severe, with a significant subsurface layer of the order of ~2 km losing its CO content very quickly before the procedure continues at a normal pace as seen in Figure 4. This change in our initial parameters greatly reduced the chances of CO survival for deep subsurfaces near the core of the nuclei, even for planetesimals we had identified as being able to maintain some CO in the other set of simulations. We note that other studies (e.g., G. Sarid &

D. Prialnik 2009) employ a pore size distribution, often with thinner tubes near the surface and wider in the interior of comets (D. Prialnik & H. Sierks 2017), which seems to be a more appropriate approach.

### 5.4.3. Implications of the Instability Timing

An important parameter of this study is the timing of the giant planet instability, which marks the beginning of our simulations. While the occurrence of the instability is widely accepted because it can explain a number of features of the solar system (see D. Nesvorný 2018 for a review), the timing remains a matter of debate. The original “Nice model” proposed an instability that was triggered relatively late (~500 Myr after the start of planet formation) to coincide with the “terminal lunar cataclysm” (R. Gomes et al. 2005; A. Morbidelli et al. 2012). Recent works have shown that the inferred cataclysm almost certainly was the result of sampling bias (P. Boehnke & T. M. Harrison 2016; N. E. B. Zellner 2017; W. K. Hartmann 2019). In addition, a late giant planet instability almost invariably would have destabilized the orbits of the terrestrial planets, assuming them to be fully formed at that time (N. A. Kaib & J. E. Chambers 2016). Various cosmochemical constraints agree that the instability must have taken place earlier than 100 million yr after the start of planet formation (A. Morbidelli et al. 2018; D. Nesvorný et al. 2018; S. J. Mojzsis et al. 2019). From a dynamical point of view, the instability could have been triggered by the dissipation of the gaseous disk (B. Liu et al. 2022), by planet-planet interactions after the dissipation of the disk (R. R. de Sousa et al. 2020), or by interactions between the giant planets and the primordial planetesimal disk (the “classic” instability trigger; K. Tsiganis et al. 2005; B. Quarles & N. Kaib 2019; R. R. de Sousa et al. 2020). A gas-disk-driven instability would have taken place on the disk dissipation timescale of a few Myr (K. E. J. Haisch et al. 2001). A self-driven instability would also have taken place quickly, within ~3–5 Myr of gas disk dispersal (R. R. de Sousa et al. 2020). A planetesimal-disk-driven instability, however, would have taken ~50 Myr to be triggered (B. Quarles & N. Kaib 2019; R. R. de Sousa et al. 2020). A gas-disk-driven or self-driven instability would have taken place while the terrestrial planets were still forming (M. S. Clement et al. 2018), whereas a planetesimal-disk-driven instability would have taken place after the main formation phase. In any case, the time between the gas disk dispersal and the triggering of the giant planet instability represents a calm period of thermal processing before the planetesimals were scattered toward their reservoirs or ejected from the solar system. As discussed above, the main difference between our results and those of C. M. Lisse et al. (2022)—who predict a lack of hypervolatiles in OC comets—comes from C. M. Lisse et al.’s (2022) assumption of a several hundred million year long stable processing phase, which is at odds with the current interpretation of early solar system evolution.

B. J. R. Davidsson (2021) examined the processing on a scenario where the instability is triggered after 15 Myr from the formation of the disk. By considering objects evolving at distances of 23 au on circular orbits, they suggested that all JFC-sized planetesimals (R = 4 km) should lose their condensed CO content in ~0.1 Myr. Their results indicate that a prolonged stay in the disk would result in significant thermal processing, pointing to the necessity of establishing a precise timeline for the instability, to get a better picture of the early processing of future
comets. More importantly, understanding the timing of planetesimal formation would allow one to more appropriately account for the effect of internal radiogenic heating, which we have neglected in this study (e.g., D. Prialnik et al. 1987; J. Haruyama et al. 1993; D. Prialnik & M. Podolak 1995; D. Prialnik 2006; U. Malamud et al. 2022).

## 6. Summary

We modeled the thermal evolution of 1518 planetesimals originating from the planetesimal disk and populating outer solar system reservoirs or escaping from the solar system to interstellar space. To this end, we used a dedicated thermal and ice evolution model coupled to the planetesimals' trajectories, obtained from two sets of N-body simulations. Our simulations showed the following.

1. Ejected planetesimals were the most processed, with KB/SD planetesimals being significantly less altered and OC planetesimals being the least processed population, on average.

2. The thermal processing during the implantation to the KB/SD and the OC or the ejection from the solar system primarily affected the pure condensed CO ice content. Less volatile species like CO₂ or amorphous H₂O ice did not present important alterations on the majority of the planetesimals.

3. Unlike previous predictions, the survival of CO ice, and therefore of other species of high volatility, is possible in a significant number of KB/SD and OC objects. This provides a natural explanation for the observed variability in the abundance ratios of CO on JFC and LPC populations and a framework for distinct drivers of cometary activity for these populations upon their return as JFCs and LPCs.

4. An intimate relation between the orbital and the thermal evolution of planetesimals was established, pointing to the necessity of taking into account the former to obtain the full picture of any past alterations on comets observed in the current epoch.

## Acknowledgments

We are grateful to Will Grundy for sharing CO data ahead of their publication. We also thank Emmanuel Lellouch and Matthew M. Knight, who reviewed this work as part of A.G.'s PhD defense and made comments that greatly improved it. 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 Pôle Scientifique de Modélisation Numérique (PSMN) of the ENS de Lyon for the computing resources. S. N.R. is grateful to the CNRS's MITI and PNP programs and also acknowledges funding from the French Programme National de Planétologie (PNP) and in the framework of the Investments for the Future program IdEx, Université de Bordeaux/RRI ORIGINS.

## Appendix

### Physical Parameters of the Thermal Evolution Model

The effective thermal conductivity accounts for reduction entailed by the porous nature of cometary material. Here we use Russell's correction factor (H. W. Russell 1935), a widely used expression in models for the thermal evolution of comets (e.g., S. Espinasse et al. 1991 and R. Orosei et al. 1999, just to name a few):

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

where f is the ratio of the solid to the radiative conductivity, which can be calculated as

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

with ϵ the material's emissivity, σ (W m⁻² K⁻⁴) the Stefan-Boltzmann constant, rₚ (m) the radius of the pores, and T (K) the temperature of the matrix. Tables 2 and 3 provide additional parameters pertaining to cometary material used in the model.
<table><caption>**Table 2**

Physical Parameters of the Nominal Object</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Units</th><th>References</th></tr></thead><tbody><tr><td>Density</td><td>ρ<sub>CO</sub></td><td>890</td><td>kg m<sup>-3</sup></td><td>R. Luna et al. (2022)</td></tr><tr><td>...</td><td>ρ<sub>CO<sub>2</sub></sub></td><td>1500</td><td>...</td><td>M. Á. Satorre et al. (2008)</td></tr><tr><td>...</td><td>ρ<sub>H<sub>2</sub>O<sub>cr</sub></sub></td><td>917</td><td>...</td><td>D. Eisenberg & W. Kauzmann (2005)</td></tr><tr><td>...</td><td>ρ<sub>H<sub>2</sub>O<sub>am</sub></sub></td><td>940</td><td>...</td><td>J. A. Ghormley & C. J. Hochanadel (1971)</td></tr><tr><td>...</td><td>ρ<sub>dust</sub></td><td>3000</td><td>...</td><td>A. Enzian et al. (1997)</td></tr><tr><td>Heat capacity</td><td>c<sub>CO</sub></td><td>35.7T - 187</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>J. O. Clayton & W. F. Giauque (1932)</td></tr><tr><td>...</td><td>c<sub>CO<sub>2</sub></sub></td><td>6.34T + 167.8</td><td>...</td><td>W. F. Giauque & C. J. Egan (1937)</td></tr><tr><td>...</td><td>c<sub>H<sub>2</sub>O</sub></td><td>7.49T + 90</td><td>...</td><td>W. F. Giauque & J. W. Stout (1936)</td></tr><tr><td>...</td><td>c<sub>dust</sub></td><td>3T</td><td>...</td><td>A. Enzian et al. (1997)</td></tr><tr><td>Conductivity</td><td>κ<sub>CO</sub></td><td>0.1</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>G. M. Muñoz Caro et al. (2016)</td></tr><tr><td>...</td><td>κ<sub>CO<sub>2</sub></sub></td><td>93.4/T</td><td>...</td><td>M. T. Mellon (1996)</td></tr><tr><td>...</td><td>κ<sub>H<sub>2</sub>O<sub>am</sub></sub></td><td>2.34 × 10<sup>-3</sup>T + 2.8 × 10<sup>-2</sup></td><td>...</td><td>J. Klinger (1980)</td></tr><tr><td>...</td><td>κ<sub>H<sub>2</sub>O<sub>cr</sub></sub></td><td>567/T</td><td>...</td><td>J. Klinger (1980)</td></tr><tr><td>...</td><td>κ<sub>Dust</sub></td><td>2.5</td><td>...</td><td>...</td></tr><tr><td>Bond albedo</td><td>A</td><td>0.04</td><td>...</td><td>...</td></tr><tr><td>IR emissivity</td><td>ε</td><td>0.96</td><td>...</td><td>...</td></tr><tr><td>Hertz factor</td><td>h</td><td>0.01</td><td>...</td><td>W. F. Huebner et al. (2006)</td></tr><tr><td>Mean pore radius</td><td>r<sub>p</sub></td><td>10<sup>-6</sup></td><td>m</td><td>...</td></tr><tr><td>Tortuosity</td><td>ξ</td><td>1</td><td>...</td><td>...</td></tr><tr><td>Local zenith angle</td><td>ζ</td><td>0</td><td>rad</td><td>...</td></tr><tr><td>Porosity</td><td>ψ</td><td>0.75</td><td>...</td><td>L. Jorda et al. (2016)</td></tr><tr><td>Nominal radius</td><td>R</td><td>5000</td><td>m</td><td>...</td></tr><tr><td>Initial temperature</td><td>T</td><td>10</td><td>K</td><td>...</td></tr></tbody></table>

<table><caption>**Table 3**

Coefficients for the Saturation Vapor Pressure of CO Given by Equation (13)</caption><thead><tr><th>Coefficient</th><th>Value</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>1.80741183 × 10<sup>1</sup></td></tr><tr><td>A<sub>1</sub></td><td>-7.69842078 × 10<sup>2</sup></td></tr><tr><td>A<sub>2</sub></td><td>-1.21487759 × 10<sup>4</sup></td></tr><tr><td>A<sub>3</sub></td><td>2.73500950 × 10<sup>5</sup></td></tr><tr><td>A<sub>4</sub></td><td>-2.90874670 × 10<sup>6</sup></td></tr><tr><td>A<sub>5</sub></td><td>1.20319418 × 10<sup>7</sup></td></tr></tbody></table>

## ORCID iDs

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

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

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

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

Nathan Kaib https://orcid.org/0000-0001-5272-5888

## References

Adachi, I., Hayashi, C., & Nakazawa, K. 1976, PThPh, 56, 1756

A'Hearn, M. F., Feaga, L. M., Keller, H. U., et al. 2012, ApJ, 758, 29

A'Hearn, M. F., Millis, R. C., Schleicher, D. O., Osip, D. J., & Birch, P. V. 1995, Icar, 118, 223

Bear, J. 1988, Dynamics of Fluids in Porous Media (New York: Dover)

Benkhoff, J., & Huebner, W. F. 1995, Icar, 114, 348

Benkhoff, J., & Spohn, T. 1991, GeoRL, 18, 261

Biver, N., Bockelée-Morvan, D., Lis, D. C., et al. 2021, A&A, 651, A25

Bockelée-Morvan, D., & Biver, N. 2017, RSPTA, 375, 20160252

Bockelée-Morvan, D., Debout, V., Erard, S., et al. 2015, A&A, 583, A6

Boehnke, P., & Harrison, T. M. 2016, PNAS, 113, 10802

Brown, G. N. J., & Ziegler, W. T. 1980, in Advances in Cryogenic Engineering, ed. K. D. Timmerhaus & H. A. Snyder (New York: Springer), 662

Chambers, J. E. 1999, MNRAS, 304, 793

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

Choukroun, M., Altwegg, K., Kührt, E., et al. 2020, SSRv, 216, 44

Clayton, J. O., & Giauque, W. F. 1932, JACChS, 54, 2610

Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icar, 311, 340

Clement, M. S., Raymond, S. N., Kaib, N. A., et al. 2021, Icar, 355, 114122

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

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 Sousa, R. R., Morbidelli, A., Raymond, S. N., et al. 2020, Icar, 339, 113605

Dones, L., Brasser, R., Kaib, N., & Rickman, H. 2015, SSRv, 197, 191

Eisenberg, D., & Kauzmann, W. 2005, The Structure and Properties of Water (Oxford: Oxford Univ. Press),

Emel'yanenko, V. V., Asher, D. J., & Bailey, M. E. 2013, EM&P, 110, 105

Enzian, A., Cabot, H., & Klinger, J. 1997, A&A, 319, 995

Espinasse, S., Klinger, J., Ritz, C., & Schmitt, B. 1991, Icar, 92, 350

Faggi, S., Lippi, M., Mumma, M. J., & Villanueva, G. L. 2023, PSJ, 4, 8

Fanale, F. P., & Salvail, J. R. 1984, Icar, 60, 476

Fanale, F. P., & Salvail, J. R. 1990, Icar, 84, 403

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

Fray, N., & Schmitt, B. 2009, P&SS, 57, 2053

Ghormley, J. A., & Hochanadel, C. J. 1971, Sci, 171, 62

Giauque, W. F., & Egan, C. J. 1937, JChPh, 5, 45

Giauque, W. F., & Stout, J. W. 1936, JACChS, 58, 1144

Gkotsinas, A., Guilbert-Lepoutre, A., & Raymond, S. N. 2023, AJ, 165, 67

Gkotsinas, A., Guilbert-Lepoutre, A., Raymond, S. N., & Nesvorny, D. 2022, ApJ, 928, 43

Golabek, G. J., & Jutzi, M. 2021, Icar, 363, 114437

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Natur, 435, 466

Gough, D. O. 1981, SoPh, 74, 21

Grundy, W., Tegler, S., Steckloff, J., et al. 2024, Icar, 40, 115767

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

Guilbert-Lepoutre, A. 2014, Icar, 231, 232

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

Haisch, K. E. J., Lada, E. A., & Lada, C. J. 2001, ApJL, 553, L153

Harrington Pinto, O., Womack, M., Fernandez, Y., & Bauer, J. 2022, PSJ, 3, 247

Hartmann, W. K. 2019, LPSC, 50, 1064

Haruyama, J., Yamamoto, T., Mizutani, H., & Greenberg, J. M. 1993, JGR, 98, 15079

Heisler, J., & Tremaine, S. 1986, Icar, 65, 13
Huebner, W. F., Benkhoff, J., Capria, M., et al. 2006, Heat and Gas Diffusion in Comet Nuclei, Vol. 133 (Noordwijk: ESA)

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

Jewitt, D. 2009, AJ, 137, 4296

Jewitt, D., Kim, Y., Mutchler, M., et al. 2021, AJ, 161, 188

Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icar, 277, 257

Kaib, N. A. 2022, SciA, 8, eabm9130

Kaib, N. A., & Chambers, J. E. 2016, MNRAS, 455, 3561

Kaib, N. A., & Volk, K. 2022, arXiv:2206.00010

Keller, H. U., Mottola, S., Skorov, Y., & Jorda, L. 2015, A&A, 579, L5

Klinger, J. 1980, Sci, 209, 271

Kossacki, K. J., Szutowicz, S. Ł., & Leliwa-Kopystyński, J. 1999, Icar, 142, 202

Kral, Q., Pringle, J. E., Guilbert-Lepoutre, A., et al. 2021, A&A, 653, L11

Lagarde, N., Decressin, T., Charbonnel, C., et al. 2012, A&A, 543, A108

Levasseur-Regourd, A.-C., Agarwal, J., Cottin, H., et al. 2018, SSRv, 214, 64

Levison, H. F., & Duncan, M. J. 1994, Icar, 108, 18

Levison, H. F., Morbidelli, A., Van Laerhoven, C., Gomes, R., & Tsiganis, K. 2008, Icar, 196, 258

Lilly, E., Jevčák, P., Schambeau, C., et al. 2024, ApJL, 960, L8

Lisse, C. M., Gladstone, G. R., Young, L. A., et al. 2022, PSJ, 3, 112

Liu, B., Raymond, S. N., & Jacobson, S. A. 2022, Natur, 604, 643

Loveless, S., Prialnik, D., & Podolak, M. 2022, ApJ, 927, 178

Luna, R., Millán, C., Domingo, M., Santonja, C., & Satorre, M. Á. 2022, ApJ, 935, 134

Malamud, U., Landeck, W. A., Bischoff, D., et al. 2022, MNRAS, 514, 3366

Malamud, U., & Prialnik, D. 2015, Icar, 246, 21

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

Meech, K. J., & Svoren, J. 2004, in Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson, AZ: Univ. Arizona Press), 317

Mellon, M. T. 1996, Icar, 124, 268

Méndez, A., & Rivera-Valentín, E. G. 2017, ApJL, 837, L1

Milani, A., Nobili, A. M., & Carpino, M. 1987, A&A, 172, 265

Mojzsis, S. J., Brasser, R., Kelly, N. M., Abramov, O., & Werner, S. C. 2019, ApJ, 881, 44

Morbidelli, A., Marchi, S., Bottke, W. F., & Kring, D. A. 2012, E&PSL, 355, 144

Morbidelli, A., Nesvorny, D., Laurenz, V., et al. 2018, Icar, 305, 262

Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes, R. 2007, AJ, 134, 1790

Muñoz Caro, G. M., Chen, Y. J., Aparicio, S., et al. 2016, A&A, 589, A19

Nesvorný, D. 2011, ApJL, 742, L22

Nesvorný, D. 2015, AJ, 150, 73

Nesvorný, D. 2018, ARA&A, 56, 137

Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117

Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, NatAs, 2, 878

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

Ootsubo, T., Kawakita, H., Hamada, S., et al. 2012, ApJ, 752, 15

Orosei, R., Capaccioni, F., Capria, M. T., et al. 1999, P&SS, 47, 839

Parhi, A., & Prialnik, D. 2023, MNRAS, 522, 2081

Prialnik, D. 2006, in IAU Symp. 229, Asteroids, Comets, Meteors, ed. D. Lazzaro, S. Ferraz-Mello, & J. A. Fernández (Cambridge: Cambridge Univ. Press), 153

Prialnik, D., Bar-Nun, A., & Podolak, M. 1987, ApJ, 319, 993

Prialnik, D., Benkhoff, J., & Podolak, M. 2004, in Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson, AZ: Univ. Arizona Press), 359

Prialnik, D., & Jewitt, D. 2022, arXiv:2209.05907

Prialnik, D., & Podolak, M. 1995, Icar, 117, 420

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

Prialnik, D., & Sierks, H. 2017, MNRAS, 469, S217

Quarles, B., & Kaib, N. 2019, AJ, 157, 67

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

Russell, H. W. 1935, Journal of the American Ceramic Society, 18, 1

Safrit, T. K., Steckloff, J. K., Bosh, A. S., et al. 2021, PSJ, 2, 14

Sarid, G., & Prialnik, D. 2009, M&PS, 44, 1905

Satorre, M. Á., Domingo, M., Millán, C., et al. 2008, P&SS, 56, 1748

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

Schorghofer, N. 2008, ApJ, 682, 697

Schorghofer, N. 2010, Icar, 208, 598

Schorghofer, N. 2016, Icar, 276, 88

Schörghofer, N., & Hsieh, H. H. 2018, JGRE, 123, 2322

Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Sci, 347, aaa1044

Simon, J. B., Blum, J., Birnstiel, T., & Nesvorný, D. 2022, arXiv:2212.04509

Steckloff, J. K., Lisse, C. M., Safrit, T. K., et al. 2021, Icar, 356, 113998

Steiner, G. 1990, A&A, 240, 533

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Natur, 435, 459

Vokrouhlický, D., Nesvorný, D., & Dones, L. 2019, AJ, 157, 181

Ward, W. R. 1974, IGR, 79, 3375

Weissman, P., Morbidelli, A., Davidsson, B., & Blum, J. 2020, SSRv, 216, 6

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

Zellner, N. E. B. 2017, OLEB, 47, 261