# On Averaging Eccentric Orbits: Implications for the Long-term Thermal Evolution of Comets

Anastasios Gkotsinas, Aurélie Guilbert-Lepoutre, Sean Raymond

► To cite this version:

Anastasios Gkotsinas, Aurélie Guilbert-Lepoutre, Sean Raymond. On Averaging Eccentric Orbits: Implications for the Long-term Thermal Evolution of Comets. The Astronomical Journal, 2023, 165 (2), pp.67. 10.3847/1538-3881/acaafd. hal-04236190

HAL Id: hal-04236190

https://cnrs.hal.science/hal-04236190v1

Submitted on 12 Oct 2023

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 CC BY 4.0 - Attribution - International License
# On Averaging Eccentric Orbits: Implications for the Long-term Thermal Evolution of Comets

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

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

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

Received 2022 October 18; revised 2022 December 8; accepted 2022 December 11; published 2023 January 24

## Abstract

One of the common approximations in long-term evolution studies of small bodies is the use of circular orbits averaging the actual eccentric ones, facilitating the coupling of processes with very different timescales, such as the orbital changes and the thermal processing. Here we test a number of averaging schemes for elliptic orbits in the context of the long-term evolution of comets, aiming to identify the one that best reproduces the elliptic orbits' heating patterns and the surface and subsurface temperature distributions. We use a simplified thermal evolution model applied on simulated comets both on elliptic and on their equivalent averaged circular orbits, in a range of orbital parameter space relevant to the inner solar system. We find that time-averaging schemes are more adequate than spatial-averaging ones. Circular orbits created by means of a time average of the equilibrium temperature approximate efficiently the subsurface temperature distributions of elliptic orbits in a large area of the orbital parameter space, rendering them a powerful tool for averaging elliptic orbits.

Unified Astronomy Thesaurus concepts: Comets (280); Comet nuclei (2160); Comet interiors (272); Computational methods (1965)

## 1. Introduction

Considering comets' long lifetimes and complex dynamical histories, studying the long-term evolution of their nuclei is an essential step in understanding their current state and activity (Gkotsinas et al. 2022). One approach to assess long-term effects is to determine their thermal histories by coupling a thermal evolution model to orbital trajectories from dynamical simulations (Raymond et al. 2020; Gkotsinas et al. 2022). This coupling is challenging both from a physical and a numerical point of view (Gkotsinas et al. 2022) and demands a certain number of assumptions, as the processes involved operate on very different timescales. Thermal evolution processes such as phase transitions, heat or gas diffusion, depending on the temperature conditions, usually take place in minutes or hours or up to a few months in the case of the heat diffusion (Prialnik et al. 2004). The dynamical evolution of comets, on the other hand, is a billion-year process that requires N-body simulations with an output frequency on the order of hundreds to thousands of years (e.g., Nesvorný et al. 2017; Sarid et al. 2019), far longer than the short timescales of thermal processes.

One of the assumptions frequently used in these types of simulations is to simplify the dynamical pathways of comets by averaging their orbits (Prialnik & Rosenberg 2009; Guilbert-Lepoutre 2012; Snodgrass et al. 2017; Gkotsinas et al. 2022). Prialnik & Rosenberg (2009) proposed to replace elliptic orbits by equivalent circular ones; the radius of which was chosen to assure the same amount of total energy over the course of an orbital period. This is reinforced in other fields, as climate modeling of putative Earth-like planets on eccentric orbits has generally found that the controlling factor determining whether a planet may retain liquid water is the total energy received

over a planet's orbit (Williams & Pollard 2002; Bolmont et al. 2016).

In this work we aim to compare and evaluate different methods for averaging elliptical orbits in the context of long-term simulations of comets' thermal evolution. In Section 2 we present the tested averaging schemes and we give a brief presentation of the thermal evolution model. In Section 3 we present the produced internal temperature distributions from the different schemes and a comparative study between them and elliptic orbits. In Section 4 we discuss the adequacy and the limits of validity of every scheme and highlight the scheme that worked better in the current context.

## 2. Methods

### 2.1. Averaging Eccentric Orbits

A variety of approaches have been proposed in order to average an elliptic orbit. The most common one is to estimate an average distance between the object and the focal point. This can be achieved in two ways: The first is to integrate the orbit equation over the true anomaly ($\theta$) throughout an orbital period to obtain a true-anomaly-averaged radius (e.g., Curtis 2014):

$$ \bar{r}_{\theta} = a\sqrt{1 - e^2}. \qquad (1) $$

The second way is to integrate the orbit equation over the average angular velocity of the object over the course of an orbital period. This approach results in a time-averaged radius, which is always smaller than the true-anomaly-averaged radius (e.g., Curtis 2014):

$$ \bar{r}_{t} = a \left( 1 + \frac{e^2}{2} \right). \qquad (2) $$

Another widely used method is to integrate over a physical parameter, such as the flux received by an object over an orbital period. This is the most physically plausible way to approach

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Subsurface temperature distributions for a layer of 10 m over a period of 1000 yr for (a) an elliptical orbit with $a = 10$ au and $e = 0.5$, (b) true-anomaly-averaged radius ($\bar{r}_{\theta}$), (c) true-anomaly-averaged flux ($\bar{r}_{\theta F}$), (d) effective thermal radius ($r_T$), (e) time-averaged flux ($\bar{r}_{tF}$), and (f) time-averaged radius ($\bar{r}_t$) equivalent circular orbits.</figcaption></figure>

the problem, as it ensures that the total energy intercepted by an object over an orbit is not modified. Using Equations (1) and (2) we can calculate a true-anomaly-averaged and a time-averaged flux (Méndez & Rivera-Valentín 2017) from which we can obtain a new set of radii, the second of which is commonly used on planetary habitability studies (e.g., Bolmont et al. 2016):

$$ \bar{r}_{\theta F} = \frac{a(1 - e^2)}{\sqrt{2 + e^2}}, \qquad (3) $$

$$ \bar{r}_{tF} = a(1 - e^2)^{1/4}. \qquad (4) $$

Recently, Méndez & Rivera-Valentín (2017) proposed a new effective thermal radius, calculated directly from the time average of the equilibrium temperature ($T_{eq}$), guaranteeing the same average equilibrium temperature over an orbital period:

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

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

where E is the complete elliptic integral of the second order.

All of the proposed expressions are simple functions of the semimajor axis (a) and the eccentricity (e). In practice, this means that for an elliptic orbit with specified orbital parameters (i.e., an (a, e) couple), the orbit-averaging technique produces a circular orbit around the focal point with an “equivalent” radius. The

differences in the averaging expressions imply different distances from the Sun. In fact for any given (a, e) couple the calculated radii are ordered as $\bar{r}_{\theta F} < \bar{r}_{\theta} < \bar{r}_{tF} < r_T < \bar{r}_t$. This means that true-anomaly-averaged flux expression ($\bar{r}_{\theta F}$) will always place an object closer to the Sun than the time-averaged radius ($\bar{r}_t$). As a consequence for the same (a, e) couple we produce different temperature profiles (panels (b) to (e) in Figure 1), raising the question “which one is better approximating the temperature distribution of the elliptic orbit?” (panel (a) in Figure 1). To answer we test these orbit-averaging schemes in the context of comets’ thermal evolution. In order to test the validity limits of each scheme in a range of orbital parameter space relevant for the inner solar system (i.e., ~3–30 au) we use a total of 110 (a, e) couples. We sample the semimajor axis range logarithmically (10^x, with x ranging between 0.5 and 1.5 with a step of 0.1) with more orbits close to the Sun, where the heating is stronger and the eccentricity is linear between 0 and 0.9 at increments of 0.1.

## 2.2. Thermal Evolution Model

We use a 1D version of the 3D thermal evolution model described in Guilbert-Lepoutre et al. (2011) to solve the heat diffusion equation in a spherical airless object:

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

where $\rho_{\text{bulk}}$ is the object’s bulk density (kg m⁻³), c is the material’s heat capacity (J kg⁻¹ K⁻¹), T is the temperature (K),
<figure><img src="image_4.png" /><figcaption>Figure 2. Temperature profiles for a subsurface layer of 10 m, ~1 Myr after the start of the simulations, for three a, e couples: (a) a = 7.94 au and e = 0.7, (b) a = 10.0 au and e = 0.5, and (c) a = 25.11 and e = 0.2. The solid black lines give the temperature profiles of elliptic orbits at perihelion, aphelion, and halfway through—time-wise—both inward and outward. The temperatures profiles for the equivalent orbits are true-anomaly-averaged radius ($\bar{r}_0$) (blue loosely dashed-dotted line), true-anomaly-averaged flux ($\bar{r}_{0F}$) (yellow dashed-dotted line), time-averaged radius ($\bar{r}_t$) (green dashed line), time-averaged flux ($\bar{r}_{tF}$) (red dotted line), and effective thermal radius ($r_T$) (purple loosely dashed line).</figcaption></figure>

$\kappa$ is the material's effective thermal conductivity (W K⁻¹ m⁻¹), and $S$ is the heat sources and sinks.

The surface boundary condition for Equation (7) is

$$ (1 - \mathcal{A}) \frac{L_{\odot}}{4\pi d_H^2} = \varepsilon\sigma T^4 + \kappa \frac{\partial T}{\partial z}, \quad (8) $$

where the received solar energy is given as a function of Bond's albedo ($\mathcal{A}$), the solar constant ($L_\odot$) in W m⁻², and the heliocentric distance ($d_H$) in au; the nucleus' thermal emission as a function of the emissivity ($\varepsilon$), the Stefan-Boltzmann constant ($\sigma$), and the temperature ($T$) in K; and the heat flux toward the interior given as a function of the surface's thermal conductivity ($\kappa$) in W K⁻¹ m⁻¹. We assume that the incident solar energy is uniformly distributed over the surface of the sphere, providing a spherical average of the energy received by the nucleus.

As the goal of this work is to compare elliptic to circular orbits, we chose a simplified setup for our model. Each comet is composed of dust without any ice, such that no phase transitions take place, removing any energy sources or sinks from Equation (7) ($S = 0$). Without them the most important parameter in the model is the effective thermal conductivity ($\kappa$), as it controls the heat diffusion toward the interior. In the current configuration it is set at $5 \times 10^{-3}$ W m⁻¹ K⁻¹, in good agreement with laboratory measurements for porous dust aggregates (Krause et al. 2011). The rest of the model's parameters are widely used averages in the published literature (Huebner et al. 2006).

It is the value of the heliocentric distance ($d_H$) that changes between the different case studies, controlling the amount of energy received at the surface of our objects. In an elliptic orbit it changes following the constant increment of the eccentric anomaly as the object moves between the apsides. In the test cases it is constant, set on the distance calculated by Equations (1)–(6).

We run a total of 660 simulations: 110 reference simulations for all the (a, e) couples with objects on elliptic orbits, serving as basis for the comparisons with the 550 simulations for objects on equivalent circular orbits created from Equations (1)–(6) for the same (a, e) couples. The simulations

run for ~1 Myr, an arbitrary period selected to allow heat diffusion inside our objects. This allows us to study the temperature differences between elliptic and circular orbits, not only at the surface but in the interior as well, and look for any accumulative or propagation effects that might be introduced during long-term simulations.

## 3. Results

Figure 1 presents the internal temperature distribution produced by the different orbit-averaging schemes for a comet on an orbit with a = 10 au and e = 0.5. For clarity only the first 1000 yr of the 1 Myr simulation are presented, but this is sufficient to notice the averaging effects. Clearly, none of the averaged orbits reproduce the heating cycle of the elliptic orbit with the subsequent passages from the perihelion to the aphelion, i.e., the seasons. Instead, as the distance from the Sun is constant, the heat diffusion is steady and uniform throughout an orbital period.

Figure 2 presents temperature profiles in the interior of our simulated comets using different orbit-averaging methods, taken near the end of our simulations (~1 Myr). The fact that the profiles are near-vertical stems from the assumption of a fixed, circular orbit with constant illumination. A profile from the reference elliptic orbit with the same orbital elements (i.e., a and e) is shown at perihelion and aphelion, and at two points halfway in between at times coinciding with a quarter of its orbital period, both on its way inward (from aphelion to perihelion) and outward (from perihelion to aphelion). The three panels examine different orbits: (a) a highly eccentric and relatively short orbit with a = 7.94 au and e = 0.7; (b) a longer, less eccentric orbit with a = 10 au and e = 0.5; and (c) a long orbit (a = 25.11 au) with low eccentricity of 0.2.

The temperature profiles from Figure 2 help us to confirm some of the previous observations:

1. The true-anomaly-averaged flux ($\bar{r}_{0F}$) better approximates the effect of the perihelion passage (even for the less eccentric orbit, where the temperature difference at the surface is ~6 K), but fails completely on the other positions both on the surface and the subsurface.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Temperature differences between the elliptic orbits at perihelion and the time-averaged radius ($\bar{r}_t$) (top row), the time-averaged flux ($\bar{\phi}_{tf}$) (middle row), and the effective thermal radius ($r_T$) equivalent orbits (bottom row) for three depths: surface (left column), 1 m (middle column), and 10 m (right column) for all the *a*, *e* couples.</figcaption></figure>

Divergences range from ~6 K in the least eccentric to ~80 K in the most eccentric orbit examined.

2. The true-anomaly-averaged radius ($\bar{r}_0$) works slightly better, with calculated surface temperatures between those of the two apsides and subsurface profiles closer to the reference ones, and with temperature differences ranging from ~20 K in the most eccentric, to almost complete convergence in the least eccentric orbits.

3. Overall, the time-averaged expressions (Equations (2), (4), and (6)) work better in all cases. Despite their failure to reproduce the high surface temperatures encountered at perihelion, especially for the most eccentric orbit, these

schemes better reproduce the cooling effect of the aphelion passage, leading to internal temperature distributions that converge to those of the reference orbit in all the examples of Figure 2.

4. The effective thermal radius ($r_T$) stands out as the scheme producing an internal temperature distribution converging almost perfectly to the reference distribution, close to the surface: below 6 m in the first two cases and 10 m in the last.

Given that the true-anomaly-averaged (or spatial) expressions failed in matching the reference simulations, we focus hereafter on the time-averaged formulas. We expand our study
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Temperature differences between the elliptic orbits at aphelion and the time-averaged radius ($\bar{r}_i$; top row), the time-averaged flux ($\bar{\phi}_{IF}$; middle row), and the effective thermal radius ($r_T$) equivalent orbits (bottom row) for three depths: surface (left column), 1 m (middle column), and 10 m (right column) for all the *a*, *e* couples.</figcaption></figure>

by examining the temperature differences from the reference elliptic orbits at the surface, 1 and 10 m below the surface at perihelion, at aphelion and halfway in between time-wise.

In Figure 3 we present the temperature differences ($\Delta T$) between the three time-averaged formulas ($\bar{r}_i$, $\bar{r}_{iF}$, and $r_T$) and the elliptic orbits at perihelion. All schemes underestimate the surface temperatures during the perihelion passage. The time-averaged radius ($\bar{r}_i$) deviates the most from the reference orbits (maximum $\Delta T$ of ~358 K or in terms of relative difference by 73%), especially for very eccentric ($e > 0.5$) and short orbits ($a < 10$ au). It is followed by the effective thermal radius ($r_T$) with a maximum divergence of ~343 K (relative difference of 70%). The time-averaged flux ($\bar{\phi}_{IF}$) presents a maximum $\Delta T$ of ~298 K (or relative difference of 61%). These relative

differences at the surface rise with the eccentricity, becoming important (~20%) above $e = 0.3$, and really significant for high eccentricities (50% for $e = 0.7$ and ~70% for $e = 0.9$) for the time-averaged and the effective thermal radius, highlighting their failure to represent the perihelion passage. The time-averaged flux ($\bar{\phi}_{IF}$) has slightly smaller relative deviations especially above $e > 0.4$ where it is constantly lower than the other time-averaged schemes by ~4%–9%.

All three schemes are more robust in the objects' interiors (middle and right panels of Figure 3). At 1 m below the surface the time-averaged radius ($\bar{r}_i$) differences converge quickly for all orbits (maximum $\Delta T$ of ~26 K) except for distant ($a > 10$ au) and highly eccentric ($e > 0.5$) orbits for which the maximum $\Delta T$ is ~49 K. The same stands for the effective
<figure><img src="image_7.png" /><figcaption>**Figure 5.** A typical Jupiter-family comet trajectory taken from Gkotsinas et al. (2022) plotted over the temperature differences between elliptic orbits and an average of the effective thermal radius ($r_T$) inward and outward for all the *a*, *e* couples. Each circle represents a point in our orbital parameter space sampling. The color code gives the scale of the temperature difference.</figcaption></figure>

thermal radius ($r_T$) only with better convergence for short ($a<10$ au) and low-eccentric orbits ($e<0.5$). For distant ($a>10$ au) and highly eccentric ($e>0.5$) orbits the problem remains but is slightly less pronounced (with a maximum $\Delta T$ of ~44 K). Long-period orbits with prolonged excursions into hot areas around perihelion allow the heatwave to advance deeper in the interior and extend the temperature differences well below the surface. In these cases underestimating the perihelion temperature remains problematic, unlike the cases of short orbits ($a<10$ au) where the heating, although more intense, takes place in a shorter period of time that is not sufficient for its diffusion in the interior (panel (a) versus panel (b) in Figure 2, for example). These differences almost completely disappear 10 m below the surface in the case of the effective thermal radius ($r_T$; $\Delta T_{max} = -1$ K for $e=0.9$; lower right panel in Figure 3) and the time-averaged radius ($\bar{r}$; with the exception of very short, $a<10$ au, and highly eccentric orbits, $e>0.7$). On the contrary, the time-averaged flux ($\bar{\Delta}T_F$) differences at the interior (middle panels in Figure 3) fail to achieve convergence in high eccentricities ($e \ge 0.5$), whether it is a short or a long orbit. In addition these deviations persist at larger depths (middle right panel of Figure 3), as there is still no convergence for highly eccentric orbits ($e \ge 0.5$) and the differences on low-eccentricity orbits are higher than those of the effective thermal radius and the time-averaged radius.

At the aphelion the temperatures are overestimated by the averaging schemes (Figure 4). The time-averaged radius ($\bar{r}_t$) works better than the other methods in the surface and the close subsurface area ($\Delta T_{max} = 20$ K versus $\Delta T_{max} = -78$ K for $\bar{r}_{tF}$ and $\Delta T_{max} = 33$ K for $r_T$ at the surface and $\Delta T_{max} = 9$ K versus $\Delta T_{max} = 50$ K and $\Delta T_{max} = 10$ K, respectively, 1 m below the surface). However, at 10 m below the surface the convergence for the effective thermal radius is almost complete ($\Delta T_{max} = -1$ K in only three short and eccentric orbits) unlike for the time-averaged radius and the time-averaged flux where there is still no convergence. For completeness, similar

representations for the temperature differences halfway through inward and outward presented in Figure 2 are given in Figures 6 and 7. As expected the temperature differences are less important comparing to those of the two apsides, but we can still observe the efficiency of the time-averaged and especially of the effective thermal radius over the time-averaged flux method.

## 4. Discussion

Overall the time-averaged schemes work better than the true-anomaly or spatial-averaged ones. This is because the true-anomaly formulas are restricted to the calculation of an average distance from the focal point. Although this seems to be a sufficient assumption, it ignores crucial information: the different time spent by an object at different distances from the focal point. In fact an object in an elliptic orbit will move much faster close to perihelion than close to aphelion, implying —in our case—more time in colder regions. This information is integrated in temporal expressions rendering them more appropriate in the approximation of elliptic orbits. With that in mind, the surface temperatures obtained from the time-averaged schemes, closer to the aphelion temperatures of an elliptic orbit, are more appropriate than an average of the perihelion–aphelion temperatures calculated by the spatial-averaged schemes.

When it comes to the time-averaged orbits, we demonstrated (Figures 3 and 4) that the time-averaged ($\bar{r}_t$) and the effective thermal radius ($r_T$) better approximate the temperature distributions of elliptic orbits with the exception of the surface temperatures at perihelion. These two distances are increasing functions of the eccentricity, whereas the time-averaged flux ($\bar{\Delta}T_F$) is a decreasing one (Méndez & Rivera-Valentín 2017). This implies that for a given (*a*, *e*) couple, the time-averaged flux will place an object closer to the perihelion, leading to a systematic overestimation of the surface temperature. On the other hand, the time-averaged and the effective thermal radius
(very close by definition; see Equations (2) and (6)) place the object closer to the aphelion accounting better for the lower temperatures reigning during the biggest part of an orbit (especially for highly eccentric ones), managing better to represent the internal temperature distribution. Interestingly, there is no clear distinction on the efficiency between these two formulas, as the time-averaged radius, as expected by its definition, works better at aphelion and halfway through both inward and outward at the surface and 1 m below, but fails to convergence as quickly as the effective thermal radius that manages to converge in all cases at maximum 10 m below the surface.

We tested the equivalent semimajor axis proposed by Prialnik & Rosenberg (2009), deriving also from a time-averaging integral ($a_c = a(1 - e^2)$). As this average distance provides the same energy per orbit as the real eccentric one, it is reliable only when orbital periods are very close. Otherwise significant deviations are observed (for very eccentric orbits for instance), with a systematic overestimation of the internal temperatures. Its validity thus remains limited to low-eccentricity orbits ($e < 0.3$) or limited timescales, as done in Guilbert-Lepoutre (2012), Snodgrass et al. (2017), and Gkotsinas et al. (2022).

When the actual orbital trajectory of a typical comet is accounted for, it is clear that orbit averaging remains a viable technique for the long-term cometary thermal evolution. Figure 5 shows the orbital changes of a simulated Jupiter-family comet during its trajectory toward the inner solar system (taken from Gkotsinas et al. 2022 with dynamical trajectories from Nesvorný et al. 2017), over an average of the inward and outward surface temperature differences for the effective thermal radius ($r_T$). With the exception of short and highly eccentric orbits ($a < 10$ au and $e > 0.5$), which are very rare, the averaging is very efficient in all the other areas of the orbital

parameter space. Even in the area of long and highly eccentric orbits ($a > 10$ au and $e > 0.5$), the temperature divergences do not overcome the ~10 K. This suggests that using the effective thermal radius to average elliptic orbits is a very effective tool with a large area of validity in the orbital parameter space and not important discrepancies outside of this area (as even for the short and highly eccentric orbits the divergence is not bigger than 20 K).

We therefore recommend the use of the effective thermal radius scheme, $r_T \approx a(1 + \frac{1}{8}e^2 + \frac{21}{512}e^4 + O(e^6))$, in long-term thermal evolution studies, as it can efficiently approximate the internal temperature distribution of an airless dusty body with low thermal inertia, such as a comet or an asteroid, evolving in eccentric orbits.

This study is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 802699). We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for computing resources. The authors would like to thank Benoit Carry for his useful comments on previous work, which provided the background for the current paper.

*Facilities: PSMN, ENS de Lyon.*

# Appendix
## Supplementary Figures

In Figures 6 and 7 we present, for completeness"the temperatures differences, as presented in Figures 3 and 4 of the main text"for the halfway points inward and outward, respectively.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Temperature differences between the elliptic orbits halfway through inward (time-wise) and the time-averaged radius ($\bar{r}_i$; top row), the time-averaged flux ($\bar{r}_{F}$; middle row), and the effective thermal radius ($r_T$) equivalent orbits (bottom row) for three depths: surface (left column), 1 m (middle column), and 10 m (right column) for all the *a*, *e* couples.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Temperature differences between the elliptic orbits halfway through outward (time-wise) and the time-averaged radius ($\bar{r}_i$; top row), the time-averaged flux ($\bar{r}_{F}$; middle row), and the effective thermal radius ($r_T$) equivalent orbits (bottom row) for three depths: surface (left column), 1 m (middle column), and 10 m (right column) for all the *a*, *e* couples.</figcaption></figure>
ORCID iDs

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

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

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

References

Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. 2016, A&A, 591, A106

Curtis, H. D. 2014, in Orbital Mechanics for Engineering Students (Third Edition), ed. H. D. Curtis (3rd ed.; Boston, MA: Butterworth-Heinemann) doi:10.1016/B978-0-08-097747-8.00002-5

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

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

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

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

Krause, M., Blum, J., Skorov, Y., & Trieloff, M. 2011, Icar, 214, 286

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

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

Prialnik, D., Benkhoff, J., & Podolak, M. 2004, in Comets II, ed. M. C. Festou & H. . Keller, Vol. 374 (Tucson, AZ: Univ. of Arizona Press)

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

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

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

Snodgrass, C., Agarwal, J., Combi, M., et al. 2017, A&ARv, 25, 5

Williams, D. M., & Pollard, D. 2002, IJAsB, 1, 61