# A Near Surface Temperature Model of Arrokoth

Orkan M. Umurhan New Horizons Science Team, Co-Investigator  SETI Institute,
189 Bernardo Way, Mountain View, CA 94043, U.S.A.  NASA Ames Research Center;
Mail Stop 245-3, Moffett Field, CA 94035, USA  Cornell Center for Astrophysics
and Planetary Sciences, Cornell University, Ithaca, NY, USA  William M. Grundy
New Horizons Science Team, Co-Investigator  Lowell Observatory, 1400 W. Mars
Hill Road, Flagstaff, AZ 86001, USA  Michael K. Bird New Horizons Science
Team, Science Affiliate  Argelander-Institut für Astronomie, Universität Bonn,
53121 Bonn, Germany.  Rheinisches Institut für Umweltforschung, Universität zu
Köln, 50931 Cologne, Germany.  Ross Beyer New Horizons Science Team, Co-
Investigator  SETI Institute, 189 Bernardo Way, Mountain View, CA 94043,
U.S.A.  NASA Ames Research Center; Mail Stop 245-3, Moffett Field, CA 94035,
USA  James T. Keane New Horizons Science Team, Science Affiliate  Jet
Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91125,
USA  Ivan R. Linscott New Horizons Science Team, Science Affiliate  Stanford
University, Stanford, CA 94305, USA  Samuel Birch Department of Earth,
Atmospheric and Planetary Science, Massachusetts Institute of Technology,
Cambridge, MA 02139  Carver Bierson New Horizons Science Team, Science
Affiliate  School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287, USA  Leslie A. Young New Horizons Science Team, Co-
Investigator  Southwest Research Institute, 1050 Walnut Street, Suite 300,
Boulder, CO 80302, USA  S. Alan Stern New Horizons Science Team, Principal
Investigator  Southwest Research Institute, 1050 Walnut Street, Suite 300,
Boulder, CO 80302, USA  Carey M. Lisse New Horizons Science Team, Science
Affiliate  Space Exploration Sector, Johns Hopkins University Applied Physics
Laboratory, 11100 Johns Hopkins Rd, Laurel, MD 20723, USA  Carly J.A. Howett
New Horizons Science Team, Co-Investigator  Southwest Research Institute, 1050
Walnut Street, Suite 300, Boulder, CO 80302, USA  Silvia Protopapa New
Horizons Science Team, Co-Investigator  Southwest Research Institute, 1050
Walnut Street, Suite 300, Boulder, CO 80302, USA  John R. Spencer  New
Horizons Science Team, Co-Investigator  Division of Space Science and
Engineering, Southwest Research Institute, Boulder, CO 80302, USA  Richard P.
Binzel New Horizons Science Team, Co-Investigator  Department of Earth,
Atmospheric and Planetary Science, Massachusetts Institute of Technology,
Cambridge, MA 02139  William B. McKinnon New Horizons Science Team, Co-
Investigator  Department of Earth and Planetary Sciences and McDonnell Center
for the Space Sciences, Washington University, St. Louis, MO 63130, USA  Tod
R. Lauer New Horizons Science Team, Co-Investigator  National Optical-
Infrared Astronomy Research Laboratory, National Science Foundation, Tucson,
AZ 85726, USA  Harold A. Weaver New Horizons Science Team, Co-Investigator
Johns Hopkins University Applied Physics Laboratory, Laurel, MD 20723, USA
Catherine B. Olkin New Horizons Science Team, Co-Investigator  Division of
Space Science and Engineering, Southwest Research Institute, Boulder, CO
80302, USA  Kelsi N. Singer New Horizons Science Team, Co-Investigator
Division of Space Science and Engineering, Southwest Research Institute,
Boulder, CO 80302, USA  Anne J. Verbiscer New Horizons Science Team, Co-
Investigator  Department of Astronomy, University of Virginia,
Charlottesville, VA 22904, USA  Alex H. Parker New Horizons Science Team, Co-
Investigator  Division of Space Science and Engineering, Southwest Research
Institute, Boulder, CO 80302, USA

###### Abstract

A near surface thermal model for Arrokoth is developed based on the recently
released $10^{5}$ facet model of the body. This thermal solution takes into
account Arrokoth’s surface re-radiation back onto itself. The solution method
exploits Arrokoth’s periodic orbital character to develop a thermal response
using a time-asymptotic solution method, which involves a Fourier transform
solution of the heat equation, an approach recently used by others. We display
detailed thermal solutions assuming that Arrokoth’s near surface material’s
thermal inertia ${\cal I}=$ 2.5 W/m-2K-1s1/2. We predict that at New Horizons’
encounter with Arrokoth its encounter hemisphere surface temperatures were
$\sim$ 57-59 K in its polar regions, 30-40 K on its equatorial zones, and
11-13 K for its winter hemisphere. Arrokoth’s orbitally averaged temperatures
are around 30-35 K in its polar regions, and closer to 40 K near its
equatorial zones. Thermal reradiation from the surrounding surface amounts to
less than 5% of the total energy budget, while the total energy ensconced into
and exhumed out Arrokoth’s interior via thermal conduction over one orbit is
about 0.5% of the total energy budget. As a generalized application of this
thermal modeling together with other KBO origins considerations, we favor the
interpretation that New Horizons’ REX instrument’s $29\pm 5$K brightness
temperature measurement is consistent with Arrokoth’s near surface material’s
being made of sub-to-few mm sized tholin-coated amorphous H2O ice grains with
1 W/m-2K-1s1/2 $<{\cal I}<$10-20 W/m-2K-1s1/2, and which are characterized by
an X-band emissivity in the range 0.9 and 1.

Classical Kuiper belt objects(250) — Small Solar System bodies(1469) — Natural
satellite surfaces(2208) — Computational Methods(1965)

##  1 Introduction

The bilobate object Arrokoth (discovered by Buie et al., 2020, as 2014 MU69)
is a flattened ($\sim$ 10km$\times$20km $\times$30km), 15.9 hour rotating,
high obliquity ($\sim 99^{\circ}$) Kuiper Belt Object (KBO) encountered by the
New Horizons spacecraft on January 1, 2019 (JD 2458485, Stern et al., 2019).
Owing to its relatively low eccentricity and inclination
($e=0.03,i=2.4^{\circ}$, respectively) and its location in the Kuiper Belt
($a=44.58$AU, and orbital frequency $\omega=6.69\times 10^{-10}$s-1), it falls
into the class of so-called cold classical KBOs (or CCKBO, for short) and, as
such, is considered a Kuiper Belt planetesimal, being perhaps one of the
oldest and relatively unprocessed relics of the solar system’s formation era
(McKinnon et al., 2020). Its surface, with low mean hemispherical albedo
($A\approx 0.063$, Stern et al., 2019; Hofgartner et al., 2021), exhibits a
relatively featureless H2O-free infrared spectrum save for the possible
presence of surface methanol ice (CH3OH) observed in absorption (Grundy et
al., 2020). Arrokoth does not appear to exhibit signs of harboring any
volatile species (Lisse et al., 2021).

Furthermore, thermal X-band emission of Arrokoth (4.2-cm wavelength, 7.2 GHz)
was observed by New Horizons’ Radio Science Experiment (REX) in both its face
on approach (low-phase angle) and look-back (high-phase angle) perspective.
Fig.1 shows a visual view of Arrokoth a few minutes prior to the look-back
phase REX scan. Analysis of the winter night side thermal emission is
consistent with a mean observed brightness temperature $T_{{\rm b,obs}}=29\pm
5$ K (Grundy et al., 2020, also Bird et al. 2022, this volume). Although New
Horizons was unable to directly measure Arrokoth’s thermal inertia, based on
other observed properties of CCKBOs it is assumed to lie somewhere between
1-10 tiu (Lellouch et al., 2013; Müller et al., 2020, also 1 tiu = 1 Wm-2K-1
s1/2), indicative of highly insulating material – at least for surface
thicknesses corresponding to diurnal skin depths ($\sim$1-2 mm). For such a
thermal inertia value applied to Arrokoth, together with an assumed low heat
capacity material, predicts a seasonal (i.e., one orbit timescale $=297.6$yr)
thermal skin depth of about 1 m (Grundy et al., 2020). Geomorphological
analysis of Arrokoth’s approach hemisphere seem to exhibit relatively bright
surface units, especially in regions close to Arrokoth’s neck (Spencer et al.,
2020). What relationship, if any, do these features have with the kinds of
surface temperatures likely on Arrokoth’s surface? Furthermore, how does
Arrokoth’s relatively bright neck region correlate to both the seasonal cycle
of received insolation and resulting surface temperatures there? Lastly, can
the REX brigthness temperature measurement be connected to and properly
reconciled with a realistic temperature model for Arrokoth?

Figure 1: View of Arrokoth (CA07) taken by the LORRI camera on its departure
trajectory approximately 10 minutes prior to the second REX scan (Grundy et
al., 2020). The lit crescent seen here would have appeared somewhat diminished
during the actual REX scan ($\sim 10^{\circ}$ shift in viewing geometry).
Figure 2: Predicted temperatures based on the simple theory of Section 2. Each
panel shows the predicted temperatures as a function of distance from Sun
($r$) and thermal inertia ($\cal I$). All models assume $A=0.06$ and
$\varepsilon=0.9$. The background flux consists entirely of the microwave
background. The panels depict (left panel) predicted summer season
temperatures $T_{{}_{s,s}}$, (middle panel) predicted interior temperatures
$T_{{}_{{\rm int}}}$, and (right panel) predicted winter season temperatures
$T_{{}_{s,w}}$. The vertical hatched lines show the location of 44 AU and
various iso-temperature lines of interest are shown with white curves.

It is thus important to generate a thermophysical temperature map of Arrokoth
to help address these questions. Stereo pair images taken by New Horizons’
Long-Range Reconnaissance Imager (LORRI) Spencer et al. (2020) have been used
to produce several cartographic data products including both a global shape
model and a relatively detailed topographic map of Arrokoth’s closest approach
day-side hemisphere. Together with this $\sim 10^{5}$ element shape model and
Arrokoth’s known orbital elements, Arrokoth’s instantaneous distance from the
Sun ($r$), as well as its subsolar latitude along every point along its orbit
are also known (Porter et al., 2018).

Although New Horizons did not attempt a direct gravity measurement of
Arrokoth, an estimate for the body’s bulk mass density may be made by applying
the gravitational slope analysis of Richardson & Bowling (2014) by using the
shape/topographic model of Spencer et al. (2020), which suggests that the
body’s mean density may be as little as $\rho=150$ kg/m3 and as high as
$\rho=650$kg/m3(McKinnon et al., 2020; Keane et al., 2020, Keane et al. 2021).
For our purposes we will adopt a nominal value of $\rho=250$ kg/m3, also see
detailed discussion in Keane et al. (2021).111There are no direct constraints
on Arrokoth’s density. The quoted values are based on an empirical trend
inferred from the observations of other small airless bodies of the solar
system. For a body of given shape and rotation, the distribution of
topographic slopes corresponds to maximum stability, i.e., the statistical
slope distribution is furthest removed from the material’s angle of repose
(Richardson & Bowling, 2014; Richardson et al., 2019). McKinnon et al. (2020)
finds this state of maximum stability by tuning the mean density to around 235
kg/m3 with the quoted 1$\sigma$ error range. We are cautioned that these are
only indirect estimates based on geophysical inference and analogy to comets
and asteroids.  Fully compacted KBOs have typical densities of about $\sim
1500$kg/m3 (e.g., see compilation found in Bierson & Nimmo, 2019), thus if
Arrokoth is made of the same stuff, then its bulk porosity $p$ reasonably
falls somewhere within the range of 0.6 to 0.85. Also, for reference we note
the only direct measured density of a comet nucleus is that for 67P/Churyumov-
Gerasimenko at around 533 kg/m3 (Pätzold et al., 2016). However we do note
that Arrokoth’s topmost layers – on the scale of a few seasonal thermal skin
depths – may have mean densities that deviate from the bulk average either due
to compaction by impacts over geologic time or due to porosity increase due to
volatile sublimation during the body’s earliest times after formation.
However, compaction appears to be less likely owing to the paucity of impact
craters on Arrokoth (Stern et al., 2019; Spencer et al., 2020; McKinnon et
al., 2020).

Another key input is the unconstrained value of Arrokoth’s surface ice’s heat
capacity at constant pressure, $C_{p}$, a quantity not directly measured by
New Horizons’ flyby of the body. For the thermal solutions developed in Grundy
et al. (2020), a value $C_{p}=350$ J/K/kg was adopted based on estimates for
cold H2O Ih ice where $C_{p}\approx 180$ J/K/kg (Castillo-Rogez et al., 2012),
and similar to adopted empirical forms like that found in Klinger (1980) and
Shulman (2004). The nearly doubled value can be rationalized on the argument
that organic impurities and/or the presence of methanol in the matrix would
substantially boost the ice complex’s heat capacity. This point is examined
further in Section 6.1.3.

With these input ingredients, together with some assumptions regarding the
surface thermal inertias, one may construct a current era temperature map of
Arrokoth’s surface, which is a vast improvement over the one presented in
Umurhan et al. (2019) and Grundy et al. (2020), which was based on a
relatively coarse $\sim 2000$ element shape model.

The main purpose of this study is detailing a general framework for generating
temperature models of Arrokoth based on this 105 facet shape model. Although
we do not directly address the questions posed earlier in this introduction –
especially connecting the predicted temperature profiles to the departure
phase REX _brightness temperature_ measurement (see companion study by Bird et
al. 2022) – the temperature solutions developed here will be employable toward
answering all of them in future studies.

The rest of this study is structured as following: Section 2 presents a simple
thermal solution in order to ground our intuition as to what to expect.
Section 3 describes the thermal solution method, which includes a discussion
of the shape model, the self-obscuration analysis of the shape model, a
detailed description of how insolation is calculated over the course of one
orbit, and finally a description of the Fourier transform based solution to
the thermal diffusion equation. Section 4 surveys our results in which we
break up our discussion into encounter day thermal properties and orbitally
averaged features. In Section 5 we develop the analysis framework that
connecting the observed brightness temperature to the thermal solutions
developed here. Section 6 discusses aspects of our results together with
several speculations including a theoretical discussion that regarding the
nature of Arrokoth’s thermal inertia, as well as a few caveats. Section 7
concludes with a short summary. The various appendices contain details of the
methods used for our thermal solutions as well as details of the radiative
transfer modeling employed here.

##  2 A Simplified Physics-Based Empirical Thermal Model

To guide our intuition we present a simple thermal model to ground
expectations. The model is one dimensional and is intended to crudely
represent a thermal response of a facet over the course of one orbit with
frequency $\omega$. Its construction is ’physics-based’ in the sense that it
is a pared-down casting of the more complete thermophysical model described in
the next section.

The model is constructed in two parts and interpreted in terms of “winter” and
“summer” time exposure given Arrokoth’s high obliquity, relatively flattened
shape, and consequent extreme polar winters and summers. Our view of this is
guided by the observation that at high summer the daylit side of Arrokoth is
close to a flat slab baking under the Sun, while at high winter the same slab
is in total darkness while at the equinoxes, Arrokoth’s surfaces are only
tangentially illuminated.

The first part describes the temperature profile resulting on the surface and
interior while the face is illuminated by the Sun with absorbed flux
$f_{{}_{{\rm eff}}}=(1-A){f}_{{}_{\odot}}$, in which ${f}_{{}_{\odot}}$ is the
peak summer solstice illumination flux, and where $A$ is the material albedo.
If $r$ is the body’s orbital radius (in AU) with $f_{{}_{{\rm SC}}}=1366$W/m2
(defined as the solar constant), then the peak solar irradiance is
${f_{{}_{\odot}}}=f_{{}_{{\rm SC}}}/r^{2}$. In this statement of the problem
$\omega=\omega_{{}_{\oplus}}r^{-3/2}$ where $\omega_{{}_{\oplus}}=1.99\times
10^{-7}$ s-1. The summer season’s subsurface temperature model,
$T_{{}_{s}}(z)$, is

$T_{{}_{s}}=T_{{}_{{\rm int}}}+\left(T_{{}_{s,s}}-T_{{}_{{\rm
int}}}\right)e^{kz},$ (1)

in which $T_{{}_{s,s}}$ and $T_{{}_{{\rm int}}}$ are the summer season surface
temperature and deep interior temperature (respectively). The thermal skin
depth scale ($\ell$) is given by $2\pi/k$ where

$k\equiv\sqrt{\rho C_{p}\omega/K},$ (2)

with $K$ being the effective material conductivity, $\rho$ the ice mass
density, $C_{p}$ the ice specific heat at constant pressure. We will often
refer to the effective thermal inertia of the medium, defined as

${\cal I}\equiv\sqrt{K\rho C_{p}}.$ (3)

The above postulated energy balance at the surface during the day entails

$\underbrace{(1-A){f}_{{}_{\odot}}+f_{{}_{{\rm bg}}}}_{{\rm received\ energy\
flux}}=\underbrace{Kk\left(T_{{}_{s,s}}-T_{{}_{{\rm int}}}\right)}_{{\rm
conducted\ thermal\ flux}}+\underbrace{\epsilon_{{\rm ir}}\sigma
T_{{}_{s,s}}^{4}}_{{\rm radiative\ losses}},$ (4)

in which $\epsilon_{{\rm ir}}$ is the material’s thermal infrared emissivity,
$\sigma$ is the Stefan-Boltzmann constant and $f_{{}_{{\rm bg}}}$ represents a
constant background illumination source.

The second part of the model describes the temperature response during the
winter season, in which the surface is illuminated only by the background
radiation field $f_{{}_{{\rm bg}}}$, which we here adopt to be just the cosmic
microwave background, i.e., $f_{{}_{{\rm bg}}}\approx f_{{}_{{\rm
cmb}}}\approx 3.15\times 10^{-6}$W/m2, corresponding to $T_{{\rm
cmb}}=2.725$K. The corresponding winter season temperature solution is
$T_{{}_{w}}(z)$ where

$T_{{}_{w}}=T_{{}_{{\rm int}}}+\left(T_{{}_{s,w}}-T_{{}_{{\rm
int}}}\right)e^{kz},$ (5)

where $T_{{}_{s,w}}$ is the night time surface temperature. The corresponding
surface energy balance, analogous to Eq. (4), is given by

$f_{{}_{{\rm bg}}}=Kk\left(T_{{}_{s,w}}-T_{{}_{{\rm
int}}}\right)+\epsilon_{{\rm ir}}\sigma T_{{}_{s,w}}^{4}.$ (6)

These coupled sets of equations are supplemented by an expression linking the
winter and summer season solutions. For this we impose the condition that all
of the energy conducted into the surface during summer returns to the surface
in winter. Therefore we have

$Kk\left(T_{{}_{s,w}}-T_{{}_{{\rm
int}}}\right)=-Kk\left(T_{{}_{s,w}}-T_{{}_{{\rm int}}}\right)$ (7)

Thus, Eqs. (4,6,7) constitute a complete set of (algebraic) equations for the
three unknowns $T_{{}_{s,w}},T_{{}_{s,s}}$ and $T_{{}_{{\rm int}}}$. From Eq.
(7) it immediately follows that the interior temperature is simply the average
of the day and night side temperatures, i.e., $T_{{}_{{\rm
int}}}=\big{(}T_{{}_{s,s}}+T_{{}_{s,w}}\big{)}/2$. The solution curves are
governed by two non-dimensional numbers. The first of these is $\eta$, given
by

$\eta\equiv\left(\frac{f_{{}_{{\rm bg}}}}{(1-A){f}_{{}_{\odot}}+f_{{}_{{\rm
bg}}}}\right)^{1/4},$ (8)

where $\eta$ is the ratio of the surface flux temperatures between the night
and day side. From this one may define averaged summertime and wintertime flux
temperatures, $\overline{T}_{s}$ and $\overline{T}_{w}$ (respectively), in
which

$\overline{T}_{s}\equiv\left(\frac{(1-A){f}_{{}_{\odot}}+f_{{}_{{\rm
bg}}}}{\epsilon_{{\rm
ir}}\sigma}\right)^{1/4},\qquad\overline{T}_{w}\equiv\left(\frac{f_{{}_{{\rm
bg}}}}{\epsilon_{{\rm ir}}\sigma}\right)^{1/4}.$ (9)

The second non-dimensional parameter $\Gamma$ is the thermal parameter
introduced by Spencer, which we hereafter refer to as the Spencer Number, and
is defined here as

$\Gamma\equiv\frac{(\rho
C_{p}K\omega)^{1/2}\overline{T}_{s}}{(1-A){f}_{{}_{\odot}}+f_{{}_{{\rm
bg}}}}\longrightarrow\frac{(\rho C_{p}K\omega)^{1/2}(\epsilon_{{\rm
ir}}\sigma)^{-1/4}}{\left[(1-A){f}_{{}_{\odot}}+f_{{}_{{\rm
bg}}}\right]^{3/4}},$ (10)

estimating the relative contribution of downward directed thermal flux and
radiative losses to space.222This non-dimensional parameter is often denoted
in the literature by the symbol “$\Theta$”, but we use $\Gamma$ here because
$\Theta$ is used elsewhere in this manuscript to denote Fourier temperature
components. The simplified non-dimensional equations that must be
simultaneously solved are

$\displaystyle
1-\small\mbox{$\frac{1}{2}$}\Gamma(\theta_{{}_{s}}-\eta\theta_{{}_{w}})-\theta_{{}_{s}}^{4}$
$\displaystyle=$ $\displaystyle 0,$ (11) $\displaystyle
1+\small\mbox{$\frac{1}{2}$}\Gamma\eta^{-4}(\theta_{{}_{s}}-\eta\theta_{{}_{s}})-\theta_{{}_{w}}^{4}$
$\displaystyle=$ $\displaystyle 0,$ (12)

where $\theta_{{}_{s}}\equiv T_{{}_{s,s}}/\overline{T}_{{}_{s}}$ and
$\theta_{{}_{w}}\equiv T_{{}_{s,w}}/\overline{T}_{{}_{w}}$. We observe that
holding all input quantities constant but only varying the orbital radius, so
long as $(1-A){f}_{{}_{\odot}}\gg f_{{}_{{\rm bg}}}$ it follows that
$\Gamma\sim r^{3/4}$, which says that the importance of surface thermal
conduction increases as an object like Arrokoth moves further away from the
Sun.

In Figure 2, we plot the predicted day and night side temperatures, as well as
the predicted deep interior temperature, for conditions representative for
those of Arrokoth, wherein we adopt $\epsilon_{{\rm ir}}=0.9$, $A=0.06$,
$\rho=250$kg/m3, and $C_{p}\approx 350$ J/K/kg. For the sake of academic
completeness, we show what the predicted values of these varied temperature
quantities for various radial positions (10 AU $<r<$150 AU) and thermal
inertias $0.1\ {\rm tiu}<{\cal I}<100$ tiu, while keeping in mind that the
likely values of Arrokoth’s ${\cal I}$ is in the range of 1-10 tiu, together
with a nominal value of $r=44$ AU. This relatively uncomplicated physics-based
empirical model predicts a range of averaged daytime surface temperatures in
the vicinity of 55-60 K, deep interior temperatures between 34 K and 38 K,
with typical averaged nighttime temperatures in the range of 10-20 K. The
winter side temperature reflects the subsurface heat flux returning to the
surface. In this respect it represents the balance between the emerging
thermal conductive flux and blackbody radiative losses of Eq. (6). We will
revisit these predictions in light of the full calculation we develop in the
next two sections.

Figure 3: Percent of sky self-obscured of the 105 merged shape model. The views shown are approximately from the approach trajectory (CA03, top panel) and departure trajectory (CA07, bottom panel). This references discussions found in both Secs. 3.1-3.2, as well as Appendix A. Table 2: Various quantities and parameters, and their adopted values where appropriate. Quantity | Value | Reference
---|---|---
Mass density of surface materials, $\rho$ | 150-650 kg/m3 | McKinnon et al. (2020); Keane et al. (2020)
Specific heat at constant pressure, $C_{p}$ | 350 J/K/kg |  Castillo-Rogez et al. (2012) \+ Grundy et al. (2020)
Thermal inertia, ${\cal I}$ | 1-10 W/m2/K/s-1/2 | Lellouch et al. (2013)
Diurnal timescale thermal skin depth, $\ell_{{\rm day}}$ |  $\sim 1$mm | Grundy et al. (2020)
Arrokoth’s semimajor axis, $r$ | 44.58 AU | Buie et al. (2020)
Orbital frequency of Arrokoth, $\omega$ |  $6.69\times 10^{-10}{\rm s}^{-1}$ (298 yr orbit) | Porter et al. (2018)
Diurnal frequency of Arrokoth, $\omega_{{\rm d}}$ |  $1.10\times 10^{-4}{\rm s}^{-1}$ (15.9 hr day) | Buie et al. (2020)
Obliquity | 99.1∘ | Spencer et al. (2020)
Mean hemispherical albedo, $A$ | $\sim 0.063\pm 0.015$ | Hofgartner et al. (2021)
thermal infrared emissivity, $\epsilon_{{\rm ir}}$ | $\sim 0.9$ | Grundy et al. (2020)

##  3 Full Body Thermal Model: Formulation

The simple equilibrium model presented in the previous section is useful for
understanding the basic thermophysics of Arrokoth, but it assumes a slab-like
surface where each facet is only in simple radiative equilibrium with the sky,
and does not see other parts of Arrokoth’s warm surface. Here we enhance the
modeling by including the detailed shape and illumination information gained
by the New Horizons flyby of Arrokoth.

Given a body’s shape model with a set of N facets labeled by “$i$” each of
whose centers are given by the vector ${\bf r}_{i}$ in the coordinate frame of
the body’s center of mass (see section 3.2 below), we solve for its
temperature profile as a function of depth normal to the surface. This is
justified on the assumption that the seasonal thermal skin depth ($\sim 1-2$m)
is small compared to the horizontal scale of each facet ($\sim 50-100$m) and
on the assumption there are no other heat sources or sinks deep in Arrokoth’s
interior after 4.56 Gyr of evolution (c.f., Lisse et al., 2021). We therefore
define the variable $\Theta_{i}(z,t)$ to be the temperature in the z-direction
underneath each facet i, where $z=0$ is the surface. We furthermore identify
$T_{i}(t)=\Theta_{i}(z=0,t)$ as the surface temperature of facet i. For this
initial study we assume the thermal conductivity to be independent of
temperature, although the methods described here are generalizable to variable
$K$. The temperature response therefore satisfies the linear heat equation

$\rho
C_{p}\partial_{t}\Theta_{{}_{i}}=\partial_{z}\big{(}K\partial_{z}\Theta_{i}\big{)},$
(13)

subject to boundary conditions, one of which is nonlinear. We assume that the
deep interior of Arrokoth no longer has any active heat sources (an assumption
likely not true during right after it was formed) and, as such, we assume the
thermal flux goes to zero with large enough depth, i.e.,

$\lim_{z\rightarrow-\infty}{K\partial_{z}\Theta_{{}_{i}}}=0.$ (14)

The nonlinearities in this problem are expressed in the upper boundary
condition, which is statement of the balance of received, emitted and interior
transmitted energies, like embodied in the similar expressions utilized in the
previous section, namely Eqs. (4,6). We do not consider the sublimation of
volatile species given no observation of any gas emission. We furthermore
assume the albedo is uniform across the surface. Thus at each facet surface we
write

$\epsilon_{{\rm ir}}\sigma T_{i}^{4}=(1-A)f_{{}_{\odot,i}}+f_{{\rm
bg}}+K\partial_{z}T_{i}+\int_{\partial S}\varepsilon\sigma
T_{j}^{4}S_{ij}d\hat{\bf s}_{j},\ \ $ (15)

where $K\partial_{z}T_{i}\equiv K\partial_{z}\Theta_{i}\big{|}_{z=0}$. We note
that the received solar insolation on facet i, $f_{{}_{\odot,i}}$, is time
dependent (see sec. 3.3). The last term in the above expression represents the
reradiated infrared radiation integrated over all facets $j$, with
infinitesimal surface element $d{\bf s}_{i}$, that are visible to facet $i$.
We assume that the infrared albedo is 0. The actual amount of reradiated
radiation received at facet $i$ is contained in the matrix $S_{ij}$, which
includes information about the relative inclination angle between each facet
pair, inverse square law effects and other mitigating factors (see further
below). We have assumed that a given facet only radiates infrared light,
contained in the frequency integrated expression $\epsilon_{{\rm ir}}\sigma
T_{j}^{4}S_{ij}$, with no contribution due to reflected light. This neglect is
justified (see Sec. 5.3.1) on account of the relatively low value of $A$. A
proper treatment of this contribution should be included in future analyses of
Arrokoth.

Eqs. (13-15) constitute the nonlinear partial differential equations to be
solved for every $T_{{}_{i}}$ and generating a solution over one orbital
period. In the following subsections we provide a high-level view of the
solution method we here employ. The methodology used here is described in part
in the Supplemental Materials of Grundy et al. (2020) and were used to derive
the solutions presented also in Umurhan et al. (2019).

Figure 4: A Lucy-Richardson stacked deconvolved image from the LORRI CA06
sequence (images LOR$\\_$0408626328 through LOR$\\_$0408626336). A single
longitude/latitide graticule centered on the spin-pole is overlaid. The
negative-spin pole faced the Sun during the New Horizons flyby. Thus, this
view corresponds to the “southern” hemisphere.

###  3.1 The Shape Model

We use the merged shape model described in Spencer et al. (2020). That model
consists of a closed three-dimensional best-fit shape derived from monoscopic
images (Porter et al., 2019), upon which a topographic surface of the spin-
pole negative face derived from stereogrammetry is fitted and merged (Beyer et
al., 2019a). The model contains 107,506 facets, which have an average area of
0.013 km2 (with a standard deviation of 0.006 km2), which corresponds to a
typical facet diameter of about 100m – this being the highest resolution
achieved during New Horizons closest approach. A flavor of the shape model is
displayed in Fig. 3. The shape model is saved as a meshfile (.obj file format)
and will be a part of a forthcoming New Horizons PDS Small Bodies Node (SBN)
release333https://pds-
smallbodies.astro.umd.edu/data_sb/missions/newhorizons/index.shtml scheduled
for late 2021 (in preparation). Fig. 4 shows a Lucy-Richardson stacked
deconvolved LORRI image of Arrokoth overlain with a single latitude-longitude
graticule. Given Arrokoth’s $>90^{\circ}$ obliquity, the approach view is from
the “southern” hemisphere, i.e., those of the negative latitudes. Because of
Arrokoth’s concave shape, a latitude-longitude coordinate representation is
impractical for our purposes as it can lead to degeneracies (see also Keane et
al, 2021). We typically refer to such conventions for qualitative discussion
purposes only. For all of our calculations this work uses a Cartesian system
centered on the body’s center of Fig. 4 (see Beyer et al., 2019b) .

Figure 5: Subsolar latitude on Arrokoth ($\lambda_{{}_{ss}}$, top panel) and
it instantaneous distance from the Sun ($r$, bottom panel) over the course of
one Arrokoth orbit and consistent with the body’s 99.1∘ obliquity. The
encounter day properties are shown with a vertical hatched line. These are
solutions based on the study found in Porter et al. (2018).  Figure 6:
Predicted surface temperatures on encounter day (January 1, 2019) on the
assumption ${\cal I}=2.54$ tiu and $\varepsilon=0.9$. The upper left shows a
view of the body on approach while the lower left depicts the view of Arrokoth
during the CA08 observation. The upper and lower right panels are the same
model but rotated to emphasize the view of the equatorial regions. Figure 7:
Histogram plots for various predicted quantities on encounter day (January 1,
2019). Model assumes ${\cal I}=2.54$ tiu and $\varepsilon=0.9$. Vertical axis
corresponds to number of facets. Of the shown panels: (far left) Distribution
of surface temperatures, (middle left) Total received flux, (middle right)
flux received as reradiation from surrounding landscape and, (far right)
surface conductive flux in which a negative value corresponds to surface
directed flow of energy from the interior. In the last of these, note the
trailing wing corresponding to the intense flux ring seen evident in Fig. 9.
Figure 8: Four views of temperature at a depth of $\sim 1$ m at encounter day.
The arrows on the colorbar show the range of temperature values at this depth.
The four views are the same orientation as those shown in Fig. 6.  Figure 9:
Four views of conductive flux ($\Lambda\equiv K\partial_{z}T_{{}_{i}}$) at
encounter day. The arrows on the colorbar show the range of flux values (also
see Fig. 7). The white line shows the transition from inward and outward
conductive flux. The four views are the same orientation as those shown in
Fig. 6. The intense surface directed ring follows the temperature terminator
(cf., Fig. 6). Figure 10: Four views of orbitally averaged surface
temperature. The four views are the same orientation as those shown in Fig. 6.
The range of temperatures are shown with arrows, see also Fig. 12. Figure 11:
Four views of absolute minimum surface temperature of each facet over the
course of one orbit. The four views are the same orientation as those shown in
Fig. 6. The range of temperatures are shown with arrows, see also Fig. 12.
Note the higher temperatures running within the equatorial zone and
occasionally outlining the putative rims of cratered morphology lying in the
tropical zone.

###  3.2 Shape Analysis: $S_{{}_{ij}}$ and ${\rm K}_{{}_{ij}}$

A shape analysis here is concerned with assessing a “who-sees-who” network,
identifying which facets $j$ of a shape are visible to a given facet $i$. This
information is then used to determine how much re-radiation is subsequently
received. We have constructed our own shape analysis algorithm, which is
detailed in Appendix A, that largely mirrors the approach and philosophy
utilized in several similarly motivated examinations of other solar system
bodies, e.g., like that used for the Moon (Gläser & Gläser, 2019), 67P (Tosi
et al., 2019), and asteroids (Rozitis & Green, 2011). The shape and transfer
model converts the continuous integral of Eq. (15) into a discrete matrix
operation, i.e.,

$\int_{\partial S}\epsilon_{{\rm ir}}\sigma T_{j}^{4}S_{{}_{ij}}d\hat{\bf
s}_{j}\rightarrow\sum_{\forall j\ {\rm visible}}\epsilon\sigma T_{j}^{4}{\rm
K}_{{}_{ij}},$ (16)

where ${\rm K}_{{}_{ij}}$ contains knowledge of the “who-sees-who” network as
well as the amount of radiative transfer communicated between element pairs.

###  3.3 Calculating Diurnal Averaged Insolation $f_{{}_{\odot,i}}$

Because the rotation rate is so short compared to the orbital time, we
calculate the diurnally averaged insolation received by each facet $i$. Figure
5 shows both Arrokoth’s subsolar latitude, $\lambda_{{}_{ss}}(t)$, and its
instantaneous distance from the Sun, $r(t)$ (in AU), over the course of one
nearly 300 year orbit (Porter et al., 2018). The instantaneous flux of solar
radiation crossing Arrokoth’s location is therefore $f_{{}_{r}}\equiv
f_{{}_{{\rm SC}}}/r^{2}$. With respect to the shape model described in the
previous section, the direction of the Sun is given by the unit vector
$\hat{\bf n}_{{}_{{\rm sun}}}$ whose individual components are

$\displaystyle n_{x,{{\rm sun}}}$ $\displaystyle=$
$\displaystyle\cos\varphi\sin\big{(}\pi/2-\lambda\big{)},$ $\displaystyle
n_{y,{{\rm sun}}}$ $\displaystyle=$
$\displaystyle\sin\varphi\sin\big{(}\pi/2-\lambda\big{)},$ $\displaystyle
n_{z,{{\rm sun}}}$ $\displaystyle=$
$\displaystyle\cos\big{(}\pi/2-\lambda\big{)},$ (17)

wherein $\varphi$ is the longitude corresponding to peak illumination (“high
noon”) on Arrokoth. Then with respect to facet $i$ with unit normal $\hat{\bf
n}_{i}$, and where high noon occurs over longitude $\varphi$ with subsolar
latitude $\lambda\rightarrow\lambda_{{}_{ss}}$, the local flux of solar
radiation received at facet is

$f_{{}_{\oplus,i}}(\varphi,\lambda)=\left\\{\begin{array}[]{cc}0,&{\rm if\
ray\ blocked\ or\ }\hat{\bf n}_{i}\cdot\hat{\bf n}_{{}_{{\rm sun}}}<0,\\\
f_{{}_{r}}\hat{\bf n}_{i}\cdot\hat{\bf n}_{{}_{{\rm sun}}},&{\rm
otherwise}.\end{array}\right.$ (18)

where the condition $\hat{\bf n}_{i}\cdot\hat{\bf n}_{{}_{{\rm sun}}}<0$ means
the Sun is below the horizon. Note that like the shape analysis described in
the previous section, we must execute ray-tracing routine to determine if
other facets block the Sun’s rays from reaching the given facet. While this
can be extremely expensive, in practice we reduce the actual number of ray-
tracing determinations by taking into account knowledge of the maximum
landscape altitude angle $\beta_{{}_{i}}$ determined in the shape analysis: if
the sun vector $\hat{\bf n}_{{}_{{\rm sun}}}$ points with an angle with
respect to the local facet horizontal that exceeds the angle $\beta_{{}_{i}}$,
then no ray-tracing calculation is needed. In other words no ray-tracing is
needed if

$\hat{\bf n}_{i}\cdot\hat{\bf n}_{{}_{{\rm
sun}}}>\cos\big{(}\pi/2-\beta_{{}_{i}}\big{)}\geq 0.$ (19)

This criterion considerably speeds up calculations. Formally then we define
$f_{{}_{\odot,i}}$ to be the daily average of $f_{{}_{\oplus,i}}$, which
involves an integral over one Arrokoth 15.9 hour day. We sample this integral
using $N_{{}_{k}}$ equally spaced longitude values $\varphi_{{}_{k}}$ ranging
from $0$ to $2\pi$. In other words we say

$\displaystyle f_{{}_{\odot,i}}(\lambda)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\hat{t}_{{\rm day}}}\int_{0}^{\hat{t}_{{\rm
day}}}f_{{}_{\oplus,i}}\Big{(}\varphi(\hat{t}),\lambda\Big{)}d\hat{t}$ (20)
$\displaystyle\approx$
$\displaystyle\frac{1}{N_{{}_{k}}}\sum_{k=1}^{N_{{}_{k}}}f_{{}_{\oplus,i}}\big{(}\varphi_{{}_{k}},\lambda\big{)}.$
(21)

In all of our calculations we choose $N_{k}=16$.

Figure 12: Histogram plots for various orbital scale quantities. Model assumes
${\cal I}=2.54$ tiu, $A=0.06$, and $\varepsilon=0.9$. Vertical axis
corresponds to number of facets. The panels are: (left panel) percent of body
self-obscuration in reference to Fig. 3, (middle left) orbitally averaged
surface temperature $T_{{}_{{\rm int},i}}$, (middle right) orbitally averaged
flux temperature $T_{{}_{{\rm a},i}}$, (right) minimum surface temperature
over the course of one orbit. Figure 13: Four views of full temperature
variation over the course of one orbit ${\Delta T}_{{}_{i}}$. The four views
are the same orientation as those shown in Fig. 6. Far upper right graph shows
a histogram of ${\Delta T}_{{}_{i}}$. Figure 14: Four views of orbitally
averaged flux temperature $T_{{}_{{\rm a},i}}$. The four views are the same
orientation as those shown in Fig. 6. The range of temperatures are shown with
arrows, see also Fig. 12.

###  3.4 Thermal Solution

We approach the solution to Eq. (13) on the assumption that the thermal
response is purely periodic, which means to say that all transients have died
away, therefore making this a so-called time asymptotic solution like used by
Titus & Cushing (2012); Schloerb et al. (2015) and White et al. (2016). We opt
for adopting this assumption and the following Fourier transform based
solution method owing to the spectral accuracy in produces, especially with
respect to the importance of preserving the total conducted energy into and
out of the interior over the course of one orbital period. As such, the
solutions are represented as a truncated Fourier series in powers of the
orbital frequency $\omega$,

$\Theta_{i}=\sum_{n=0}^{N}\Theta_{i,n}e^{k_{n}z}e^{-i\omega nt}+{\rm c.c.},$
(22)

in which

$k_{n}=\left(\frac{i\rho C_{p}\omega n}{K}\right)^{1/2},\qquad{\rm
Re}(k_{n})>0,$ (23)

where $N$ is an integer that is about the sampling rate of the received flux
$f_{{}_{\odot,i}}$ (see further below). For each value of $n$, the functional
forms in Eq. (22) are the exact solutions of the Fourier Transform of Eq. (13)
where the amplitudes $\Theta_{i,n}$ are unknown and must be determined from
enforcing boundary conditions. Eq. (22) automatically solves the lower thermal
boundary condition Eq. (14). We observe that $\Theta_{i,0}$ is half of the
deep interior temperature ($\equiv T_{{\rm int},i}$), and so long as the
thermal forcing is not aperiodic we note that $\equiv T_{{\rm int},i}$ is time
independent for a given facet, but may vary from facet to facet. The advantage
of this approach is that we do not have to directly calculate subsurface
temperatures numerically as this information is encoded in the vertical
structure function $e^{k_{n}z}$ for each Fourier mode $n$. For the sake of
academic completeness we describe the rudiments of the solution approach in
Appendix B, and discuss further in the Appendix C how this approach may be
generalizable to configurations in which the conductivity changes value with
depth, as might be the case if one is confronted with a physical situation in
which very low thermal inertia refractories sit atop an otherwise highly
conductive solid ice layer below (see also section 5). This multi-layer model
was in fact used in generating the Arrokoth thermal solutions based on the
1962 facet model found on the Small Bodies Node PDS data repository
(https://pds-
smallbodies.astro.umd.edu/data_sb/missions/newhorizons/index.shtml), and is
the reason we have included detailing this kind of layered solution approach
in Appendix C.

With Eq. (22) in hand, it is straightforward to assess the surface thermal
flux ($\Lambda_{{}_{i}}$)

$\Lambda_{{}_{i}}\equiv
K\partial_{z}T_{{}_{i}}=K\sum_{n=1}^{N}k_{n}\Theta_{i,n}e^{-i\omega nt}+{\rm
c.c.}$ (24)

The final stage involves determining the coefficients $\Theta_{i,n}$, which,
once inverse Fourier tranformed, gives the time series of $T_{i}$ over the
course of one cycle $\omega$. The coefficients are determined using an
iterative procedure sketched in the following: for notational ease we identify
$f_{{}_{i}}\equiv\epsilon_{{\rm ir}}\sigma T_{{}_{i}}^{4}$, then the boundary
condition is

$f_{{}_{i}}=(1-A)f_{{}_{\odot,i}}+K\partial_{z}T_{{}_{i}}+{\rm{K}}_{ij}f_{{}_{j}}.$
(25)

The zeroth order solution is the solution to the above with the conductive
term set to zero. Thus

$f_{{}_{i}}^{(0)}=\Big{(}\mathtt{I}_{{}_{ij}}-{\rm
K}_{{}_{ij}}\Big{)}^{-1}(1-A)f_{{}_{\odot,i}},$ (26)

which involves a matrix inversion – with $\mathtt{I}_{{}_{ij}}$ the identity
matrix – followed by determining the first iterate temperature solution
$T_{{}_{i}}^{(1)}$ to

$f_{{}_{i}}^{(1)}=\epsilon_{{\rm
ir}}\sigma\left(T_{{}_{i}}^{(1)}\right)^{4}=f_{{}_{i}}^{(0)}+K\partial_{z}T_{{}_{i}}^{(1)}.$
(27)

In this equation $f_{{}_{i}}^{(0)}$ is known at all points along the orbit,
thus it becomes a matter of solving for $T_{{}_{i}}^{(1)}$. This is achieved
by performing a fast Fourier Transform of Eq. (27) into frequency domain
followed by using a Newton-Raphson routine 444 For example,
https://en.wikipedia.org/wiki/Newtons_method. to solve for the Fourier
components $\Theta_{i,n}^{(1)}$. The received flux, as well as the solution of
the surface temperature field, is well represented with a $N=30$ Fourier mode
decomposition. We have assessed the robustness of these $N=30$ results against
solutions generated up to $N=300$ finding characteristic differences in the
temperature field (for any given facet) differing by less than 1 part in 105.
For a given face, we demand for convergence that successive changes in each
$\Theta_{i,n}^{(1)}$ be less than $10^{-4}$, which is usually achieved in 6 or
less sub-iterations. All higher iterations $m\geq 2$ involve calculating

$\displaystyle\epsilon_{{\rm ir}}\sigma\left(T_{{}_{i}}^{(m)}\right)^{4}$
$\displaystyle=$ $\displaystyle
K\partial_{z}T_{{}_{i}}^{(m)}+(1-A)f_{{}_{\odot,i}}+{\rm
K}_{{}_{ij}}f_{{}_{j}}^{(m-1)},\ \ \ $ (28) $\displaystyle f_{{}_{i}}^{(m)}$
$\displaystyle=$ $\displaystyle\epsilon_{{\rm
ir}}\sigma\left(T_{{}_{i}}^{(m)}\right)^{4},$ (29)

where, as before, the first of the above equations is solved in Fourier
frequency space, while the second takes that solution and expresses it the
time domain, calculating the matrix product ${\rm
K}_{{}_{ij}}f_{{}_{j}}^{(m)}$ in the time domain, and then feeding the result
forward to the next iteration, $m\rightarrow m+1$. Note that while this
procedure involves one large matrix multiplication per iteration, it is not
the main bottle-neck in the calculation, but rather the Newton-Raphson stage,
for which values of $m\geq 2$ tends to converge in 1-2 sub-iterations owing to
the closeness of the initial guess to a converged solution. For this
particular system, in which the characteristic Spencer number $\Gamma$ is low,
we find that satisfactory convergence is achieved in one ($m=1$) iteration.

This calculation is run in parallel using four 3.1 GHz Intel Cores (i7). The
temperature solution for each facet (including all sub-iterations) takes less
than 0.004 cpu seconds/processor for $N=30$ and 0.09 seconds/processor for
$N=300$. We find that it takes approximately 40 minutes to generate a full
solution with one iteration for the $10^{5}$ facet model when $N=300$, whereas
it takes a little over 3 minutes when $N=30$. All results displayed in this
study are done with $N=300$.

We close this section by defining the orbitally averaged _flux temperature_
$T_{{}_{{\rm a},i}}$, which is, for each facet, the orbital average of the
solution to Eq. (15) in the infinitely insulating limit (i.e., $K\rightarrow
0$). Given the preceding discussion, this amounts to the solution to

$\epsilon_{{\rm ir}}\sigma T_{{}_{{\rm
a},i}}^{4}=\frac{1}{P}\int_{0}^{P}f_{{}_{i}}^{(0)}\Big{(}\lambda(t)\big{)}dt,$
(30)

in which $P$ is Arrokoth’s orbital period, and where $f_{{}_{i}}^{(0)}$ –
given in Eq. (26) – is understood to be a time-dependent function of the
subsolar latitude $\lambda(t)$. This quantity is a proxy for the total
received insolation for each facet. Given Arrokoth’s high obliquity, we expect
$T_{{}_{{\rm a},i}}$ to be relatively higher in the regions experiencing
extreme polar winters and summers ($|\lambda|\gtrapprox 10^{\circ}$) compared
to Arrokoth’s perpetually diurnal equatorial zones ($|\lambda|\lessapprox
10^{\circ}$). Henceforth, we designate as _polar region_ or _polar zone_
latitudes that satisfy $|\lambda|\gtrapprox 10^{\circ}$, while ascribing
_tropical zone_ or _equatorial zone_ to latitudes satisfying
$|\lambda|\lessapprox 10^{\circ}$ (e.g., as in the usage of Earle et al.,
2017).

##  4 Full Body Thermal Model: Results

###  4.1 Encounter Day Views

In Fig. 6 we present an example temperature map of Arrokoth viewed on
encounter day (January 1, 2019 or Julian Date 2458485) on the assumption that
${\cal I}\approx 2.5$ tiu, which, in turn, corresponds to an effective Spencer
Number $\Gamma_{{}_{{\rm eff}}}\approx 0.0082$ (with
$\overline{T}_{{}_{s}}\approx 59$K), where we follow the definition found in
Eqs. (9-10) and adopting
$\overline{f}_{{}_{\odot}}=f_{{}_{\odot}}\big{(}r=44{\rm AU}\big{)}$. On this
day, the subsolar latitude ($\lambda_{{}_{ss}}$) was
$\lambda_{{}_{ss}}\approx-61.86^{\circ}$ and Arrokoth was $\approx 43.24$ AU
from the sun (see also Fig. 5). Four views are shown in Fig. 6: The upper left
panel shows the approach view, where the simulated body specific subspacecraft
coordinates are at $\varphi=-46.9^{\circ},\ \lambda=76^{\circ}$, while the
lower left panel shows the departure view with nominal subspacecraft
coordinates $\varphi=46.9^{\circ},\ \lambda=284^{\circ}$. The latter view
position closely approximates the view from the spacecraft when the REX (CA08)
scan was made.

For these parameters the predicted surface temperatures ranges from 10 K to 60
K. The distribution of temperatures, in terms of facet counts, is shown in the
first panel of Fig. 7. A cursory view of the upper left panel of Fig. 6
displays the propensity for high temperatures (in the vicinity of 57-60K) over
most of the flattened and exposed real estate of the encounter hemisphere
ranging from the south pole down to $\lambda=-10^{\circ}$ or so. The view from
the REX observation point (lower left panel of Fig. 6) shows the surface
temperature is dominated by the winter side values (around 12-13 K) with a
sliver of high temperatures along the top outer rim. The temperature
transitions from the winter side lows to the summer side highs through the
equatorial region, with the transition occurring across a narrow band of
tropical estate.

Fig. 8 shows a predicted temperature map at approximately 1m below Arrokoth’
surface. The input parameters are as those shown in Fig. 6. The temperatures
at this depth range from about 40K down to 25K. Notably, the temperatures at
this depth are slightly cooler at high summer side latitudes than they are at
latitudes approaching the equator (the differences being 3-4 K). This is to be
expected since tropical zone latitudes never go into polar night over the
course of one orbit, which means that the temperature there never get nearly
as cold as they do in the poles (see further below). The winter side
temperatures, on the other hand, do not show this slight inversion. The reason
for this is likely due to the winter-side part of the shape model not being
nearly as flattened as the imaged encounter hemisphere – but we caution
against drawing strong conclusions here owing to the true uncertainty of the
topography of the unseen hemisphere.

Fig.9 shows the instantaneous daily averaged conductive flux at the surface at
the day of the encounter. Once again the view and input parameters are the
same as those shown for the two prior figures. The figure confirms one’s
intuition, in which the thermal flux is interior directed on the lit summer
side of Arrokoth ($\sim$ 0.002 W/m2, i.e., ), while the conductive flux is
surface directed on Arrokoth’s winter side ($\sim$ -0.0015 W/m2). We have
delineated where the transition occurs from outward to inward flux as well.
One notable feature is the appearance of a relatively high intensity ring of
surface directed thermal flux ($\sim$ -0.003 W/m2) along the temperature
terminator. This ring coincides with the location of the pre-sunrise edge of
Arrokoth’s surface, where the temperature gradient along the surface is
greatest (e.g., compare against the surface temperatures shown in Fig. 6). We
find this ring to closely follow the illumination terminator over the course
of Arrokoth’s solar revolution.

Fig. 7 shows several facet count histograms for various quantities at
encounter day. These include a distribution of surface temperatures to be
compared with Fig. 6 and surface fluxes of Fig. 9. The intense ring of surface
fluxes appears as a long tail below the distribution mode at
$\Lambda\approx-0.0015$W/m2. We have also shown the distribution of instaneous
total received illumation flux at every facet, which includes both direct
received sunlight as well as illumination received from reradiation from other
“seen” facets. The distribution shows a strong mode at around 0.6 W/m2, while
including a large number of facets receiving nearly zero flux (i.e., only
CMB). We have separated out from the total received flux only the flux
received from surface reradiation (labelled as flux received from surrounding
landscape). The irradiation received in the form of surface re-radiation
($\sim 0.02$ W/m2) amounts to less than 5% of the total irradiation. Once
again this appears to be due to Arrokoth’s relatively muted local relief
despite its otherwise flattened global shape. These trends are consistent with
those based on low-order 1962 facet model reported in Grundy et al. (2020).

As a final reflection, the behavior here ought to be compared to temperature
predictions done for a body with a neck but much more marked topography (e.g.,
67P/Churyumov-Gerasimenko), would be interesting in its own right, in terms of
the role of re-radiation, shadowing, etc. (Hu et al., 2017).

###  4.2 Orbital Timescale Features

Adopting ${\cal I}=2.54$ tiu, as in the prior section, we display various
orbital timescale properties. In Fig. 10 we showcase the orbitally averaged
surface temperature of each facet. This quantity also corresponds to the
asymptotic interior temperature (i.e., $T_{{}_{{\rm int},i}}$) descriptive of
subsurface positions exceeding several thermal skin depths. For all regions of
the body except for a narrow equatorial band ($\sim\pm 10^{\circ}$), the deep
interior temperature is in the low 30K-35K range. On the other hand,
$T_{{}_{{\rm int},i}}$ in the equatorial band is closer to 40K. These trends
make sense given that polar winters last long enough for all the received and
thermalized summer time solar energy to radiate away. Since the equatorial
“tropical” zone always sustains a diurnal insolation pattern, the received
solar insolation over those regions never sufficiently radiates away to lower
the surface temperatures nearly as much before the Sun rises again. This
effect is most dramatically illustrated in Fig.11 showing the absolute minimum
temperature any one facet experiences over the course of one orbit. Almost the
entirety of Arrokoth’s surface, including substantial parts of the tropical
zone, experiences minimum temperatures in the range of 10K-15K. But, a very
thin sliver of real estate following the equator and occasionally outlining
cratered morphology regions found within the tropical zone, show minimum
temperatures as high as 30K. The distribution of these temperatures are shown
in Fig. 12.

The interior depths of the equatorial zone are thus warmer than the polar
region despite the fact that these parts receive overall less insolation over
the course of one orbit. This is clear when viewing the predicted flux
temperature as defined in Eq. (30). $T_{{}_{{\rm a},i}}$ is a proxy for the
total irradiation received on any given facet and is shown in Fig. 14. As an
orbital average, the tropical zones receive far less solar insolation, with
$T_{{}_{{\rm a},i}}\sim 40$K, than the polar regions, where $T_{{}_{{\rm
a},i}}\sim 45$K (see also third panel of Fig. 12). Yet despite this, the
averaged surface temperature in the tropics is much higher than the poles.

##  5 Connecting to REX observed brightness temperatures

Figure 15: Predicted brightness temperatures based on subsurface radiative
transfer analysis for $\big{<}E_{\rm eff}\big{>}=1$. Several values of ${\cal
I}$ shown. White contours denote several values of electric skin depth. Note
the weak dependence on the $\varepsilon^{\prime}$ (donoted on graphs as
Re$(\varepsilon)$). Figure 16: Predicted brightness temperatures $T_{b}$ for
the CA08 observation. For input values of $\varepsilon^{\prime}\equiv{\rm
Re}(\varepsilon)$ and $\big{<}E_{{\rm eff}}\big{>}$ predicted $T_{b}$ shown as
a function of $\varepsilon^{\prime\prime}\equiv{\rm Im}(\varepsilon)$ and
${\cal I}$: (left column) $\varepsilon^{\prime}=1.04$, (right column)
$\varepsilon^{\prime}=3.50$. Successive rows corresponding to decreasing
values of $\big{<}E_{{\rm eff}}\big{>}$. Black contours denote $T_{b}=29$K
while dotted black contours represent $T_{b}=24$K and 34K levels. The right
column might be considered representative of pure methanol ice of zero
porosity, while the left column could be the same with porosity of 60%.

The REX radio flux density measurements made during the CA08 observing
sequence, together with the shape model for Arrokoth, leads to an Arrokoth-
disk averaged x-band brightness temperature $T_{b,{\rm obs}}=29\pm 5K$ (Grundy
et al., 2020). This result has been reaffirmed in a more comprehensive
analysis (Bird et al. 2022). The brightness temperature depends on various
properties of Arrokoth’s near surface materials. In order to calculate a model
prediction for the brightness temperature that we here simply call $T_{b}$, it
is necessary to perform a radiation transfer analysis that relates Arrokoth’s
interior kinetic temperature $T(z)$ – as detailed in previous sections – and
its X-band thermal and refractive properties to $T_{b,obs}$. In other words,

$T_{b,{\rm obs}}=T_{b}\Big{(}{\cal
I},\varepsilon^{\prime},\varepsilon^{\prime\prime},\big{<}E_{{\rm
eff}}\big{>}\Big{)},$ (31)

in which $T_{b}$ depends on knowledge of the material’s thermal inertia, the
value of the real and imaginary part of its X-band dielectric
constant/permittivity $\varepsilon=\varepsilon\prime+i\varepsilon\prime\prime$
555Note, hereafter all values of $\varepsilon$ are scaled in units of vacuum
permittivity $\varepsilon_{0}\approx 8.85\times 10^{-12}$F/m, where F is in
units of Farads., and its effective X-band emissivity $\big{<}E_{{\rm
eff}}\big{>}$. These 4 a priori unknown parameters may themselves have some
kind of depth-dependence, but for our purposes here we take them to be
constants. In principle, there exists a relationship between ${\cal
I},\varepsilon^{\prime},\varepsilon^{\prime\prime}$, and $\big{<}E_{{\rm
eff}}\big{>}$ that produces values of $T_{b}$ equal to $T_{b,{\rm obs}}$.
Anticipating our discussion in section 6.1, at best we will be able to only
circumscribe this relationship between the four unknown parameters.

Utilizing the approach detailed in Section 3, we develop a suite of near
surface thermal solutions for several values of the thermal inertia, i.e., in
the range $0.5\ {\rm tiu}<{\cal I}<275\ {\rm tiu}$ on 35 equally spaced values
along a logarithmic scale. Together with our adopted values of $\rho$ and
$C_{p}$, this range of ${\cal I}$ values corresponds to a characteristic
orbital timescale skin depth values,

$\ell_{\rm orb}\equiv{\cal I}\Big{/}\rho C_{p}\sqrt{\omega},$ (32)

that fall into the range $22\ {\rm cm}\lessapprox\ell_{\rm orb}\lessapprox
120\ {\rm m},$ where $\ell_{\rm orb}$’s definition is based on the setting
$n=1$ into the absolute value of the inverse of the thermal wavenumber defined
Eq. (23). It is sufficient to use $N=30$ Fourier modes to generate this
solution array (see section 3.4). We develop solutions in which for depths
below 10 meters from the surface the thermal conductivity transitions to a
fixed constant value for porous H2O ice (e.g., $K=1$W/m/K, Klinger, 1980,
corresponding to a value ${\cal I}=295$ tiu – also see sec. 3.4 and Appendix
C).

With these temperature solutions in hand, we then employ the radiative
transfer solution method outlined in de Kleer et al. (2021) to produce a set
of predicted $\tilde{I}_{\nu}$ that finally relate to $T_{b}$ via Eq. (31).
For the transfer solutions we consider the following ranges of the remaining
parameters: $1\leq\varepsilon^{\prime}\leq 3.5$ at 70 equally spaced values on
a linear scale, $10^{-4}\leq\varepsilon^{\prime\prime}\leq 10$ at 75 equally
spaced values on a logarithmic scale, and four values for, $\big{<}E_{{\rm
eff}}\big{>}=0.7,0.8,0.9,1.0$. These parameter choices follow those also
considered in Bird et al. (2022). Further details of the solution method as
implemented are found in Appendix D. We note that the electrical skin depth,
$\delta_{{\rm elec}}$, approximately relates to $\varepsilon$ via

$\delta_{{\rm
elec}}=\frac{\lambda}{2\pi}\frac{\sqrt{\varepsilon^{\prime}}}{\varepsilon^{\prime\prime}},$
(33)

in the $\varepsilon^{\prime\prime}\ll 1$ limit of its formal definition found
in Eq. (D4). Finally, we do not account for Fresnel reflection at the
interface.

In our companion study, Bird et al. (2022, this issue) also develop $T_{b,{\rm
obs}}$ predictions based on the thermal solutions. The differences lie in the
approaches taken to solve the radiative transfer problem: As each spacecraft
visible facet $i$ has its own unique vertical temperature profile $T_{i}(z)$,
this work develops a facet-by-facet solution to the transfer problem followed
by taking a facet area weighted average across the visible disk of Arrokoth to
derive $T_{b,{\rm obs}}$. On the other hand, Bird et al. (2022) develop a
simpler estimate for $T_{b,{\rm obs}}$ based on a single effective vertical
temperature profile $\overline{T}(z)$, which is constructed as a facet-area
weighted average of all $T_{i}(z)$, and a surface averaged facet norm with
associated cosine of the spacecraft inclination angle ${\cos\theta_{t}}$.
Formally speaking the two approaches are not necessarily commensurate – as the
average of products is not necessarily equal to the product of averages (also
see discussion in Appendix D – but in this case here we find that the results
happen to be in mutual agreement.

We observe several trends. Most prominent is that for given values of ${\cal
I}$ there is only a weak dependence of on $\varepsilon^{\prime\prime}$ as
$\varepsilon^{\prime}$ varies from 1 to 3.5 (see Fig. 15). Fig. 16 displays
contour levels of several predicted $T_{b}$ as a function of
$\varepsilon^{\prime\prime}$ and ${\cal I}$ for several fixed values of
$\varepsilon^{\prime}$ and $\big{<}E_{{\rm eff}}\big{>}$. Except for
$\big{<}E_{{\rm eff}}\big{>}=0.7$, the $T_{b}=29$K contour exists for all
input parameter values considered. In section 6.1 we further discuss and
interpret the implications of these solutions, especially with respect to
known candidate materials. However we can make one clear observation here that
solutions in which $\big{<}E_{{\rm eff}}\big{>}=0.7$ appear to be ruled out as
the $T_{b}=29$K contour fails to even register under these circumstances and,
at best, only the low end of the REX observation is predicted ($\sim 24$K) but
does so at very low (nearly transparent) values of
$\varepsilon^{\prime\prime}(<10^{-3})$ for ${\cal I}\lessapprox 40$tiu, while
the predicted $T_{b}=24$K contour is plausible only for ${\cal I}>100$tiu at
considerably larger upper bound values for $\varepsilon^{\prime\prime}$, i.e.,
in the range of 0.1 and 1.

We pose another question: Given $\varepsilon^{\prime},{\cal I}$ and
$\big{<}E_{{\rm eff}}\big{>}$, what value of $\varepsilon^{\prime\prime}$
yields a predicted $T_{b}$ during the CA08 observation? We define this
critical value as $\varepsilon_{c}^{\prime\prime}$ and we derive these values
using interpolation methods based on the solutions determined above. Fig. 17
displays $\log\varepsilon_{c}^{\prime\prime}$ as a contour plot for
$\big{<}E_{{\rm eff}}\big{>}=0.9,1.0$. We have left out doing this exercise
for the $\big{<}E_{{\rm eff}}\big{>}=0.7$ and $0.8$ cases because $T_{b}=29$K
registers as a possibility only for a limited range of ${\cal I}$ values for
the latter case and none at all for the former case (see Fig. 16). It is
notable that in both cases $\varepsilon_{c}^{\prime\prime}<0.05$ for values of
${\cal I}<70$ tiu. This analysis also predicts that the material is very
nearly transparent in the X-band (i.e., $\varepsilon_{c}^{\prime\prime}\approx
0.001$ or less) for values of $2\ {\rm tiu}\ \lessapprox{\cal I}\lessapprox
40$tiu. We revisit this in the next section.

Figure 17: Value of $\varepsilon_{c}$” that yields $T_{b}=29$K during CA08 for
$\big{<}E_{{\rm eff}}\big{>}=0.9,1.0$. For these values of $\big{<}E_{{\rm
eff}}\big{>}$ this analysis predicts that the material is very nearly
transparent in the X-band ($\varepsilon_{c}^{\prime\prime}\lessapprox 0.002$)
for $1\ {\rm tiu}\ <{\cal I}<40$tiu. Figure 18: The disk averaged “deep”
temperature $T_{{\rm deep}}$ as a function of ${\cal I}$, predicted for CA08.

A question that will be addressed later is if $T_{b,{{\rm obs}}}$ might be
explained by X-band radiation emanating from well below the orbital time scale
skin-depth. To help answer this question, we define $T_{{\rm deep}}$ to be the
visible-facet weighted observed X-band brightness temperature at
$\big{<}E_{{\rm eff}}\big{>}=1$ for values of the dielectric skin depth –
$\delta_{{\rm elec}}$ defined in Eq. (D4) – greatly exceeding $\ell_{{\rm
orb}}$. In other words, it corresponds to a visible-facet weighted disk
averaged deep subsurface temperature sampled by nearly transparent X-band
radiation sampling the deep interior. 666We are careful to note that there is
a limit here as well. One might imagine the X-band radiation is so transparent
through Arrokoth’s interior that it could also sample its sunlight side.
However this would suggest unrealistically low values of
$\varepsilon^{\prime\prime}\ll 10^{-6}$. From this interpretation it follows
that

$T_{b}\Big{(}\big{<}E_{{\rm eff}}\big{>},\delta_{{\rm elec}}\gg\ell_{{\rm
orb}}\Big{)}\rightarrow\big{<}E_{{\rm eff}}\big{>}\cdot T_{{\rm deep}}.$ (34)

In Fig. 18 we plot $T_{{\rm deep}}$ as a function of ${\cal I}$. Inspection of
the figure clearly shows that values of $\big{<}E_{{\rm eff}}\big{>}=0.7$ are
clearly precluded from predicting even $T_{b}=24$K on the low ${\cal I}$ end,
while possibly permitting $T_{b}=29$K on the high ${\cal I}$ end (i.e., $>150$
tiu). We return to discussing this in the next section.

##  6 Discussion

###  6.1 On the thermophysical properties of Arrokoth’s subsurface materials

There are no independent measures of the four input parameters that go into
these thermal solutions constructed here. Unfortunately, there is very little
one can conclude with certainty about their values or properties based on the
New Horizons’ single REX scan (see also Bird et al. 2022). Partly given its
laboratory data availability and partly given that planetesimals are thought
to be substantially composed of H2O ice, we mainly focus on the property of
water ices despite there being no direct evidence for it on the surface. This
H2O ice which may also have a tholin covering. We save for future analysis the
possibility for other less-well laboratory studied materials to explain the
REX observation.

The primary concern here is that there is a relationship between
$\varepsilon^{\prime},\varepsilon^{\prime\prime},{\cal I}$, and
$\big{<}E_{{\rm eff}}\big{>}$ that yield a predicted $T_{b}=29K$ corresponding
to the brightness temperature measured during the CA08 REX observation. One
example relationship is the quantity
$\varepsilon^{\prime\prime}=\varepsilon_{c}^{\prime\prime}\left(\varepsilon^{\prime},{\cal
I},\big{<}E_{{\rm eff}}\big{>}\right)$ developed in the previous section and,
for example, displayed in Fig. (17). Aside from Arrokoth harboring surface
methanol, there is no other information about the body’s surface or subsurface
composition. As such, we are here relegated to speculating about the possible
values of these quantities in relation to other known bodies of the outer
solar system. In this section we consider various facets of the problem and we
conclude by offering what we consider to be our favored interpretation.

####  6.1.1 Permittivity/Dielectric Constant

Perhaps most confounding for this analysis are the uncertain values of
$\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$, the latter of which
strongly controls the material’s attenuation of the X-band signal through the
subsurface medium especially when $\varepsilon^{\prime\prime}\ll 1$, e.g., see
Eq. (D4). We consider several possibilities in order of increasing
speculation:

  * •

Methanol, likely a significant constituent of Arrokoth’s near surface, has a
value of $\varepsilon^{\prime}\approx 3.52\pm 0.05$, but with no corresponding
measurement of $\varepsilon^{\prime\prime}$ (Le Gall et al., 2016). This same
study also reported a value of $\varepsilon^{\prime}\approx 3.42\pm 0.05$ for
the similar hydrocarbon ethanol. One can estimate an effective
$\varepsilon^{\prime}$ on the assumption a certain ice fraction of vacuum
porosity ($p$ hereafter) following the Maxwell Garnet empirical formula for
vacuum mixtures discussed in Bird et al. (2022). If, for example, the upper
layer is composed of pure methanol with $p=0.6$, then it would have an
effective $\varepsilon^{\prime}\approx 1.04$ corresponding to the left column
of solutions shown in Fig. 16, while the $p=0$ solution would be the same
figure’s right column.

  * •

Based on laboratory work of Paillou et al. (2008), which examined the
dielectric properties of various tholin-like materials, Bird et al. (2022)
compiled a range of possible effective $\varepsilon$ values of several types
of porous tholins based on the same prescription for mixtures. For porosity
$p$ in the 0.60 to 0.80 range: powdered tholins exhibits values of
$\varepsilon^{\prime}$ from $1.05$ to $1.01$, with
$\varepsilon^{\prime\prime}$ values ranging from $1\times 10^{-3}$ to $5\times
10^{-4}$; while for two types of compact tholins examined
$\varepsilon^{\prime}$ correspondingly ranges from $1.5$ to $1.2$ together
with $\varepsilon^{\prime\prime}$ falling somewhere between $4.2\times
10^{-3}$ and $1.5\times 10^{-3}$.

  * •

Arrokoth’s subsurface materials might resemble the “dirty” ice thought to
comprise the nucleus of comets like 67P. Heggy et al. (2012) have compiled
several known measurements of $\varepsilon=\varepsilon_{{\rm mix}}$ of matrix
ice with variable dust-to-H2O mixtures with mixture fraction $0<\phi_{{\rm
mix}}<1$.The measurements were made at wavelengths a factor 10-50 larger than
the X-band. From these they compile an empirical relationship for
$\varepsilon_{{\rm mix}}^{\prime}$ and $\varepsilon_{{\rm
mix}}^{\prime\prime}$ as a function of $\phi_{{\rm mix}}$, as well as a weak
dependence on $T$, For $T=30$K and $p=0.5$ application of their relationship
predicts for $\phi_{{\rm mix}}\ll 1$ that $\varepsilon_{{\rm
mix}}^{\prime}\approx 1.42$ and $\varepsilon_{{\rm mix}}^{\prime\prime}\approx
10^{-4}$, while for $\phi_{{\rm mix}}\rightarrow 1$ $\varepsilon_{{\rm
mix}}^{\prime}\approx 2.90$ and $\varepsilon_{{\rm mix}}^{\prime\prime}\approx
4\times 10^{-2}$. These values will be somewhat smaller in bulk for higher
porosities at fixed $\phi_{{\rm mix}}$. Similarly, Brouet et al. (2015)
compiled $\varepsilon^{\prime}$ values for an ice mixture comprised of H2O and
the lunar regolith simulant JSC-1A. For porosities $p\approx 0.6$ they find
$1.5<\varepsilon^{\prime}<2.25$, depending on mixture fraction (e.g., see
Figure 10 of Brouet et al., 2015); however, no values for
$\varepsilon^{\prime\prime}$ were reported.

####  6.1.2 X-band Emissivity

Bird et al. (2022) compile a known list of the radio emissivity of the Kronian
moons at $\lambda=2.2$cm based on Cassini RADAR observations. The color and
albedos of these objects vary widely, with corresponding emissivity range of
roughly $0.6<\big{<}E_{{\rm eff}}\big{>}<1$. However, if one focuses on the
darkest of these objects (i.e., those with the lowest albedos comparable to
Arrokoth) we find, for example, Iapetus’ leading dark side has $\big{<}E_{{\rm
eff}}\big{>}\approx 0.87$ (Le Gall et al., 2014) while for Phoebe
$\big{<}E_{{\rm eff}}\big{>}=0.92$ (Ostro et al., 2006). If these bodies’s
surface materials are analogs to Arrokoth’s, then supposing $\big{<}E_{{\rm
eff}}\big{>}=0.9$ for it would not be unreasonable.

In reference to our preceeding discussion, it seems $\big{<}E_{{\rm
eff}}\big{>}=0.7$ is highly unlikely. One might argue by inspecting the deep
temperature solutions shown in Fig. 18 that $T_{b}=29$K is achievable for
$\big{<}E_{{\rm eff}}\big{>}=0.7$ on the high ${\cal I}$ end (i.e., the
product $\big{<}E_{{\rm eff}}\big{>}\cdot T_{{\rm deep}}$ yielding 29 K), but
cross-referencing that against the $\big{<}E_{{\rm eff}}\big{>}=0.7$ row of
Fig. 16 would imply that this is feasible only for very large electric skin-
depths, corresponding to $\varepsilon^{\prime\prime}\ll 10^{-4}$. Given the
discussion in the previous section we know of no plausible materials that
would exhibit such highly transparent qualities. Performing the same cross-
comparison a similar line of reasoning applies for rejecting $\big{<}E_{{\rm
eff}}\big{>}=0.8$ solutions for values of ${\cal I}<20$ tiu. However we see
from the penultimate row of Fig. 16 that $T_{b}=29$K is achievable for
$\big{<}E_{{\rm eff}}\big{>}=0.8$ as it corresponds to physically reasonable
values of $10^{-4}<\varepsilon^{\prime\prime}<10^{-3}$ in the thermal inertia
range $20\ {\rm tiu}<{\cal I}\lessapprox 300$ tiu. Based on the discussion, of
the following section we also lean toward treating this case as unlikely (but
not ruled out) on account of identifying plausible materials with such high
bulk thermal inertias at such low temperatures and porosities.

####  6.1.3 Thermal Inertia

TNOs and KBOs as analogs. The 1-10 tiu ${\cal I}$ range of values reported in
Lellouch et al. (2013) for KBOs corresponds to daytime measurements that probe
their surface materials down to their diurnal skin depths, which are in the
vicinity of a few mm. Such low thermal inertias, which ought roughly be
considered body averages, can be interpreted as indicating the presence of a
very porous matrix of weakly loaded grains whose grain-grain radii of contact
($a$) are much smaller than the grain size (radius $R_{g}$). Since the ratio
$h\equiv a/R_{g}$ (sometimes known as the Hertz factor) controls grain-to-
grain heat conduction, small values of $h$ – perhaps due to weak loading –
coupled to high porosity conditions leads to thwarted effective heat
conduction (see further below, and also Howett et al., 2010; Ferrari & Lucas,
2016; Ferrari, 2018). Therefore, it might not be unreasonable to suppose such
fluffy low thermal inertia material extends much deeper than a few mm on
small-sized KBOs like Arrokoth as little material compaction is expected under
such low gravities ($\sim 1$mm/s2) over its natural history – unless its
surface experienced sufficient impacting events, but this appears to be ruled
out on the basis of the relatively few observed impact features (Stern et al.,
2019; Spencer et al., 2020; McKinnon et al., 2020).

Iapetus and Phoebe as analogs. While daytime observations were used to derive
thermal inertias in the previous section, more direct methods have been used
for both Iapetus and Phoebe who, at least on their surface, could be
considered analogs for Arrokoth: The three bodies share similar Bond albedos
($\sim 0.05$ Ostro et al., 2006) and Phoebe, which likely sources Iapetus’
dark side materials, may originate from the same population of KBOs as
Arrokoth (Johnson & Lunine, 2005; Castillo-Rogez et al., 2012) 777Although
recent challenges to this view have been raised based on the similarities of
Phoebe’s spectral properties to the C-type class of asteroids (Hartmann, 1987;
Matson et al., 2009), as well as open questions based on internal thermal
evolution modeling of its early history (Castillo-Rogez et al., 2019).. In
particular, Iapetus’ leading side’s very dark top soil (especially that on
Cassini Regio, CR hereafter) has been argued to be sourced from the Phoebe
ring, which is itself sourced from Phoebe via micrometeorite impact gardening
of its surface (Verbiscer et al., 2009). Thus considering Iapetus and Phoebe’s
surface materials together as an analog for Arrokoth is not wholly
unjustified, at the very least.

Based on observations made by Cassini’s Composite Infrared Spectrometer
(CIRS), estimates for the thermal inertias of several of the Kronian
satellites have been published (e.g., Flasar et al., 2005; Howett et al.,
2010; Rivera-Valentin et al., 2011). Current thermal inertia estimates for
Iapetus’s darkened side based on CIRS data are $6\ {\rm tiu}\lessapprox{\cal
I}\lessapprox 21\ {\rm tiu}$ (Howett et al., 2010), with a more recent
narrower range at $11.0\ {\rm tiu}\lessapprox{\cal I}\lessapprox 14.8\ {\rm
tiu}$ (Rivera-Valentin et al., 2011). For comparison, Iapetus’ lighter
trailing side surface materials is modeled from CIRS data to have $12\ {\rm
tiu}\lessapprox{\cal I}\lessapprox 33\ {\rm tiu}$ (Howett et al., 2010), and
possibly more tightly constrained at $15\ {\rm tiu}\lessapprox{\cal
I}\lessapprox 25\ {\rm tiu}$ (Rivera-Valentin et al., 2011). Analysis of CIRS
data provides estimates for Phoebe’s surface thermal inertias around ${\cal
I}=30$tiu (Spencer et al., 2004; Flasar et al., 2005).

It is also important to note that based on thermal modeling leveraged against
Cassini RADAR data, which was used to observe CR’s subsurface at 2.2 cm, Le
Gall et al. (2014), report higher thermal inertias of at least around ${\cal
I}\approx 50$ tiu and possibly even ${\cal I}>200$ tiu. Given that the
electric skin-depths in the radio are vastly larger than in the infrared, this
higher ${\cal I}$ estimate is likely characteristic of surface materials far
deeper than 20cm. Given the likelihood that this deeper (presumably) H2O ice
is both compacted and in cubic crystalline form, we do not consider it as a
plausible analog for Arrokoth (see further below).

The thermal inertia of amorphous H2O ice grains. Even though New Horizons did
not directly detect H2O ice on Arrokoth’s surface, theoretical considerations
including global disk evolution modeling indicate that planetesimals formed
out beyond the water snowline should be made up of significant amounts either
H2O ice grains or more complex silicate grains with H2O locked in within them
(e.g., Estrada et al., 2016, and many others). Perhaps for Arrokoth these H2O
ice grains are encased by a thin cover of methanol or other hydrocarbons.
Moreover, a recent study of the carbonaceous chondrites NWA5717 and Allende
(Simon et al., 2018) reveal that these primitive bodies are comprised of cm-
scale aggregates of sub-mm chondrules (with average diameters of the
distribution mode in the range 0.15-0.70 mm). These authors find that a
variety of subgroups of particles (characterized by different compositions
and/or lithologies) manifest similar size distributions and aggregate
character, which lead them to posit that the same size sorting process might
be operating across the entire early planetesimal formation phase of the solar
nebula.

These considerations taken together lead to the not-unreasonable proposal that
small KBO bodies like Arrokoth are composed of fluffy aggregates of H2O
-_dominated_ sub-mm particles. Assuming H2O to be the dominant component of
individual grains and assuming that the grains were assembled into Arrokoth in
situ, then it is also not unreasonable to propose that they are and have
always been in an amorphous crystalline form. This is particularly compelling
since the amorphous to cubic crystalline transition in type 1h H2O ice becomes
important on the timescale of the solar system’s existence once temperatures
exceed $70$K (Efimov et al., 2011) – a condition unlikely to have been
possible at 40-45 AU where CCKBOs reside. For reference we note that amorphous
H2O ice has a thermal conductivity that is lower than its cubic-phase
counterpart by an order of magnitude or more when compared under the same
thermophysical conditions (Andersson & Suga, 1994).

Motivated by a need to explain Mimas’ estimated thermal inertia Ferrari &
Lucas (2016) have suggested that the low ${\cal I}$ observed of a large array
outer-solar system bodies (Centaurs, TNOs, etc.) might be rationalized as
arising from primitive highly porous grains of amorphous H2O ice grains. We
entertain this possibility for Arrokoth starting with a brief theoretical
primer on the matter following Ferrari & Lucas (2016):

On an airless body, the effective thermal inertia of a fluffy aggregate of
icy-particles of size $R_{g}$ and intrinsic ice-grain density $\rho_{{\rm
ice}}$ is posited in the form

${\cal I}=\sqrt{(1-p)\rho_{{\rm ice}}C_{p}(T)(K_{{\rm r}}+K_{{\rm c}})},$ (35)

The effective conductivity, $K_{{\rm eff}}\equiv K_{{\rm r}}+K_{{\rm c}}$, is
the sum of the grain-to-grain heat transport via radiation ($K_{{\rm r}}$) and
the thermal conduction from grain to grain across shared contact points
($K_{{\rm c}}$). For $K_{{\rm r}}$ there are several known approaches (e.g.,
Shoshany et al., 2002; Piqueux & Christensen, 2009) but here we adopt the
physically tractable model form of Gundlach & Blum (2012) where

$K_{{\rm r}}=8\epsilon\sigma
T^{3}R_{g}\phi_{r}(p),\qquad\phi_{r}(p)\equiv\frac{e_{1}p}{1-p},$ (36)

where $\epsilon$ is the IR emissivity. Grossly speaking, the expression
represents the transport of blackbody radiative energy across the chasm
separating grains of similar size (i.e., $\sim R_{g}$). The porosity dependent
radiation exchange factor $\phi_{r}(p)$ accounts for the added effective
distance this radiation must travel due to increasingly rarefied interstitial
structure of low porosity grain aggregates. $e_{1}=1.34$ is a structural
constant of the model.

The contact mitigated conductivity can be described with
Johnson–Kendall–Roberts theory of contact mechanics (JKR hereafter, Johnson et
al., 1971),

$K_{{\rm C}}=h\phi_{c}(p)K_{{\rm A}},$ (37)

where $\phi_{c}(p)$ is another porosity dependent factor taking into account
reduced conduction due to the increase of vacuum space as $p$ approaches 1.
Here we assume the form $\phi_{c}(p)=2(1-p)$, which is an averaged
approximation based on empirical data across a range of porosities as reported
in Gusarov et al. (2003) (also see Ferrari & Lucas, 2016). $K_{{\rm A}}$ is
the intrinsic conductivity of an amorphous H2O ice grain that can be estimated
on theoretical grounds to be given by

$K_{{\rm A}}=C_{p}\rho_{{\rm ice}}v_{{\rm ph}}\lambda_{\rm{mol}}/4,$ (38)

where $v_{{\rm ph}}$ is its phonon speed, $\lambda_{{\rm mol}}$ is the typical
molecule separation in the crystal. Adopting this theoretical estimate checks
out well against experimental data at temperatures slightly higher than that
of interest here ($\sim 70$ K, Andersson & Suga, 1994). Absent experimental
verification in the colder temperature range relevant here, we assume Eq. (38)
is valid for our purposes. We treat H2O ice’s heat capacity – which is assumed
to be the same in its cubic crystalline and amorphous forms – according to its
known temperature dependence, for which Shulman (2004) has developed an
empirical form for based on experimental data in the range of 0 and 273K (see
Eq. 4 of Shulman, 2004). The Shulman approximation is an improvement over the
often used one attributed to Klinger (1980), in which
$C_{p}(T)=7.49(T/1{\rm{K}})+90$ J/m/kg.

The grain-grain contact radius $a$ reflects the elastic deformation resulting
from adhesive “pinching” forces related to a grain’s intrinsic surface tension
plus the deformation that occurs under a given load. For the kinds of low
gravities of concern here, it is expected that adhesive forces will dominate
overburden load down to at least a few meters on Arrokoth. Thus for amorphous
water ice grains under negligible load the Hertz factor simplifies to

$h=\frac{a}{R_{g}}=\left(\frac{L_{{\rm ec}}}{R_{g}}\right)^{1/3},\qquad
L_{{\rm ec}}\equiv\frac{9\pi(1-\nu_{p}^{2})\gamma}{4E},$ (39)

where $E$ is the ice grain’s Young’s modulus, $\nu_{p}$ is the Poisson ratio,
and $\gamma$ is the ice grain’s surface tension. The elastocapillary length
scale $L_{{\rm ec}}$ nominally denotes deformation scale at which restorative
elastic forces balance adhesive ones arising from surface tension (e.g.,
Jagota et al., 2012; Style et al., 2013).

Fig. 19 shows theoretical predictions of ${\cal I}$ for several values of
porosity in which $\epsilon=\epsilon_{{\rm IR}}=0.9$ with additional parameter
inputs found in Table 3. For putative grain sizes $0.1$mm$<R_{g}<1$mm and
temperature range $20$K $<T<50$K, the thermal inertia of amorphous H2O ice
with $p=0.6$ broadly falls into the range $5$ tiu $<{\cal I}<25$ tiu, while
for $p=0.8$ ${\cal I}$ it lies within $2$ tiu $<{\cal I}<12$ tiu.

These figures also show that for a given $T$ and $p$ there are values of
$R_{g}$ where ${\cal I}$ is minimized. This is easy to understand as the
radiative conductivity $K_{\rm r}$ has a linear dependence on $R_{g}$ while
the contact limited ice conductivity $K_{\rm C}$ diminishes as $R_{g}^{-1/3}$:
the former becomes more efficient as radiation travels across longer distances
in the matrix as $R_{g}$ increases while for fixed $L_{{\rm ec}}$ the latter
becomes more effective because smaller grain sizes means that grain-grain
contact area approaches the surface area of the grains themselves (i.e.,
$h\rightarrow 1$). This means that for given inputs ($T,p,$ etc.) that ${\cal
I}\geq{\cal I}_{m}$, where

${\cal I}_{m}=1.325\left[(1-p)C_{p}\rho_{{\rm ice}}\Big{(}8\epsilon\sigma
T^{3}K_{{\rm A}}^{3}L_{{\rm ec}}\phi_{{\rm r}}\phi_{{\rm
c}}^{3}\Big{)}^{1/4}\right]^{1/2},$ (40)

at corresponding grain size $R_{g}=R_{m}$

$R_{{\rm m}}=\left(\frac{L_{{\rm ec}}^{1/3}\phi_{{\rm c}}K_{{\rm
A}}}{24\epsilon\sigma\phi_{{\rm r}}T^{3}}\right)^{3/4}.$ (41)

The numerical factor on ${\cal I}_{m}$ comes from the minimization procedure
and is an approximation to $\sqrt{3^{1/4}+3^{-3/4}}$. As can be gleaned from
Fig. 19 we note that these extreme values occur for grain sizes of interest
based on origins considerations (sub-to-few mm size) within the range of
plausible ice porosities relevant for Arrkoth ($0.6<p<0.8$). Fig. 20 shows
${\cal I}_{m}(T,p)$, and we observe that the for the relevant temperature and
porosity ranges, ${\cal I}\leq{\cal I}_{m}\approx$1-5 tiu are not possible,
thus forming a reasonable lower bound for ${\cal I}$. We also observe from the
form of $R_{{\rm m}}$ in Eq. (41), as well as from scanning the solutions
shown in Fig. 19, that ${\cal I}_{m}$ achieves its minimum in the 20-50 K
temperature range for grain sizes between $0.7$mm$<R_{g}<$ 5mm.

Table 3: Various thermophysical and related quantities for amorphous H2O ice.

Quantity | Value | Note
---|---|---
$E$ |  $\sim 2.4\times 10^{10}$ Pa | near $T=50$K
$\nu_{p}$ | $0.33$ |
$\gamma$ |  $0.07$-$0.37$ J/m2 | range of measurments
$\rho_{{\rm ice}}$ |  $940$ kg/m3 |
$v_{{\rm ph}}$ |  $2500$ m/s |
$\lambda_{{mol}}$ |  $5\times 10^{-10}$ m |

For source references see Ferrari & Lucas (2016).

Figure 19: Predicted thermal inertia of amorphous ice as a function of grain
size and temperature for two porosities: (top) $p=0.6$, (bottom) $p=0.8$ and
$\gamma=0.37$J/m2 and $\epsilon=\epsilon_{{\rm IR}}=0.9$. White contours
highlight key ${\cal I}$ values for reference  Figure 20: Minimum thermal
inertia for amorphous H2O ice as a function of $T$ and $p$ with
$\epsilon=0.9$. Under given conditions this quantity implies that ${\cal
I}\geq{\cal I}_{m}(T,p)$. Note also this quantity is extremized at around
$p\approx 0.32$.

####  6.1.4 A speculative synthesis

Despite the analysis of the previous sections, we can at best offer likely
possible candidates for Arrokoth’s near surface materials that are consistent
with the $T_{\rm b,obs}=29\pm 5$K measurement made by the REX instrument. This
speculation is proffered on the basis of certain uncontroversially sober
assumptions that: (1) Arrokoth’s surface materials and/or layers are
uncompacted retaining a porosity in the range $0.6<p<0.8$ consistent with
Arrokoth’s mean density $\sim 250-350$ kg/m3 – on the assumption that its
constituent grains are typical of the mean densities of relaxed KBOs ($\sim
1600$ kg/m3, e.g., Castillo-Rogez et al., 2012; Bierson & Nimmo, 2019); (2) As
being a member of the CCKBOs, Arrokoth’s surface temperature probably never
exceeded 70 K at any time during its natural history owing to its stable
location in the Kuiper Belt; (3) the X-band value of the imaginary part of the
dielectric permittivity of Arrokoth’s subsurface materials is no less than
that of unpolluted H2O ice, i.e., $\varepsilon^{\prime\prime}>10^{-4}$.

By analogy to Phoebe and Iapetus’ leading side’s dark materials – mainly owing
to their shared low albedos – we favor the assumption that Arrokoth’s X-band
emissivity is in the range $\big{<}E_{{\rm eff}}\big{>}=0.9-1$. We rule out
$\big{<}E_{{\rm eff}}\big{>}$ values of 0.7 on the basis that it would
correspond to ice with $\varepsilon^{\prime\prime}$ well below $10^{-4}$. We
deem $\big{<}E_{{\rm eff}}\big{>}=0.8$ unlikely for similar reasons on the low
${\cal I}$ end ($<20$ tiu). $\big{<}E_{{\rm eff}}\big{>}=0.8$ seems unlikely
on the high ${\cal I}$ end ($>80$ tiu) despite its being a possible solution
(e.g., see Fig. 16) as it would correspond to ${\cal I}$ for either very low
or very high porosity H2O ice (both amorphous or cubic). According to Fig. 19
the high porosity ice grain solution that predicts ${\cal I}>50$ tiu would
correspond to boulder sized solid ice grains, which is difficult to produce if
the grains are formed during the era of planetesimal formation (e.g., Estrada
et al., 2016). The low porosity ice grain solution that yields such high
thermal inertias would correspond to sub-micron sized amorphous H2O ice grains
that, while possible, is hard to reconcile with the meteoritic record and
planetesimal formation models.

If Arrokoth’s subsurface materials are composed mainly of mm-scale amorphous
H2O ice grains, then its thermal inertia is likely below 10-20 tiu since
higher thermal inertias are not reasonable to expect for amorphous H2O ice at
the porosities of interest here. Of course, thermal inertias in the range
50-100 tiu are possible, but only if the H2O ice is in cubic form, which is a
possibility only if the near surface materials have experienced sustained
periods with temperatures higher than 70-80 K (see extensive discussion in
Ferrari & Lucas, 2016). As mentioned above, high values of ${\cal I}$ are
possible for amorphous ice if the grains are like boulders in size where
thermal transport is governed by radiation, but such large constituent surface
material is hard to justify from an origins perspective. The same theoretical
concerns suggest that Arrokoth’s thermal inertias are no less than $2-5$ tiu,
thus forming a reasonable lower bound for ${\cal I}$ as lower values are not
predicted for amorphous H2O ice with porosities $0.6<p<0.8$.

Taken together we see from inspecting the low $\varepsilon^{\prime}$ column of
Fig. 16 corresponding to its first two rows ($\big{<}E_{{\rm
eff}}\big{>}=0.9,1.0$) that for $1\ {\rm tiu}<{\cal I}<10$ tiu, the $T_{\rm
b,obs}=29\pm 5$K prediction occurs in the range
$10^{-4}<\varepsilon^{\prime}\lessapprox 10^{-2}$, which is consistent with
permittivity values for H2O ice in the GHz range as briefly surveyed in
section 6.1.1 and discussed further in Bird et al. (2022). Amorphous H2O ice
of mm-scale is therefore a strong candidate material to account for the
Arrokoth’s subsurface materials.

However, the same solution survey cannot rule out tholins – hydrocarbons
produced from plasma discharging in mixtures of methane and nitrogen – as
candidate materials as they are expected to be ubiquitous reddening agents of
the outer solar system (Cruikshank et al., 2005a). Consideration of laboratory
data of tholins summarized in section 6.1.1 as well as in Bird et al. (2022)
indicates that high porosity tholins exhibit
$10^{-4}<\varepsilon^{\prime}\lessapprox 10^{-3}$, which is also consistent
with $T_{{\rm b,obs}}$ for both $\big{<}E_{{\rm eff}}\big{>}=0.9,1.0$.
Unfortunately, thermal inertias for tholins are not well constrained and while
the $T_{{\rm b,obs}}$ measurement is consistent with the permitted range of
thermal inertias for amorphous H2O ice (i.e., 1-20 tiu) for $\big{<}E_{{\rm
eff}}\big{>}=1.0$, the full range of thermal inertias are permitted for
$T_{{\rm b,obs}}$ in the case where $\big{<}E_{{\rm eff}}\big{>}=0.9$ (second
row of Fig. 16). We suspect that the surface ice is likely an admixture of
amorphous H2O ice and tholins, with the latter forming a small percentage of
the composite, much in line with suggestions made for the make-up of the
Kronian satellites (Cruikshank et al., 2005b). If the dominant constitutent is
the H2O component, then we suspect the lower range of permissible thermal
inertias characterizes the $\big{<}E_{{\rm eff}}\big{>}=0.9$ case too.

We observe that in the recent study by Ferrari et al. (2021), where they
examined the thermal response of the Kronian moons – but not including Iapetus
and Phoebe – that a nominal upper bound conductivity of $K_{e0}\sim
0.001$W/m/K was adopted for their model tholin-covered amorphous H2O ice
grains. If such a tholin-H2O ice complex were present on Arrokoth, it would
correspond to an upper bound value $\sqrt{(1-p)\rho_{{\rm ice}}K_{e0}}\approx
11$ tiu.

###  6.2 On the comparison between simple theory and full results

In Section 2 we constructed a simplified algebraic theory for the predicted
temperatures on the winter and summer sides. It is of interest to see how well
this theory works against our predictions based on the more detailed
sophisticated model we have developed here. Examining Fig. 5 shows that at
encounter day Arrokoth was in its late southern hemisphere spring phase of its
orbit. Given that it was near the southern summer solstice, it is useful to
compare the distribution of the predicted encounter day surface temperatures
against the simplified theory of Sec. 2. In Fig. 21 we re-plot the surface
temperature distribution shown in the first panel of Fig. 7, where we have
distinguished between surface facets in the polar and tropical zones. Surface
facets are designated as belonging to the polar zone if there exists at least
one Arrokoth day in which the sun never rises, while tropical zones are those
facets in which the sun rises daily. The distribution shows two modes. The low
end mode corresponds to approximately 12.5K, which is the predicted value for
$T_{{}_{s,w}}\approx 12.2$K based on the results shown in Fig. 2 for ${\cal
I}\approx 2.54$ tiu (the low end of likely ${\cal I}$ values as discussed in
the previous section). The high temperature mode shows that the peak
temperature sits at around 58.5K. The simple theory predicts a typical summer
time temperature of $T_{{}_{s,s}}\approx 60$K, which is only slightly higher
than high end mode seen in the distribution, and where this discrepancy is
likely due to comparing the simple theory against a solution for a time
slightly shy of solstice (cf., Fig. 5). We conclude that the simple theory of
Sec. 2 is overall useful especially with respect to predicting the
temperatures on both the summer and winter sides near around solstice.

Figure 21: A more detailed view of encounter day surface temperatures.
Comparisons between simple theory predictions are shown: vertical line at
$T\sim 60$K corresponds to the “summer” side predictions of the simple theory
while the vertical line at $T\sim 12.5$K corresponds to the “winter” side
predictions of the same.

###  6.3 On the value the time-asymptotic solution method

The traditional approach to solving Eqs. (13-15), while similarly involves
precomputing the shape model analysis and insolation profile, calls for
developing an initial value solution to Eq. (13). Accurately forward time
evolving the heat equation can be treacherous and often time-consuming,
especially if the parameters of the problem make the system numerically stiff
(requiring short time-stepping), although there are many recent methods
developed in the astrophysical literature that alleviate this computational
overhead while preserving accuracy and stability in various problems involving
diffusion operators (e.g., Meyer et al., 2012).

The main issue, however, is developing a solution to a periodically driven
problem that is “time-asymptotic”, in the sense that the long-time state of
the system – certainly after tens of millions of orbits since formation
(McKinnon et al., 2020) – is also periodic with the external driver’s
periodicity (if, however, out of phase, etc.). Arrokoth is likely in such a
time-asymptotic state given that its orbital elements are not thought to have
experienced strong chaotic or aperiodic episodes, at least not in recent
history (Porter et al., 2018). Thus, the problem we have here lies in
initiating the initial value problem with the right initial data that sets the
solution off into an exactly periodic transient-free state. But this is
something which, in general, is not known _a priori_. In practice, such
initial value problems are setup up with some initial profile (satisfying
boundary conditions), and subsequently forward evolved until the time-periodic
asymptotic state is eventually reached. But this may take time to achieve
(sometimes dozens of forcing periods), which can be expensive. The Fourier
transform method utilized here, as well as in Titus & Cushing (2012); White et
al. (2016), while slightly costly due to the fast Fourier transform itself,
immediately admits an exactly periodic solution, with no transient structure.
There is no forward time stepping involved and the resulting solutions are
accurate to spectral precision, which makes this solution method attractive
from the standpoint of generating reliable solutions.

###  6.4 On thermally driven mechanical erosion, a speculation

Beyond circumstantial evidence for methanol based on spectral fitting to LEISA
data obtained by New Horizons (Stern et al., 2019; Grundy et al., 2020), what
are Arrokoth’s near surface materials and what are their thermo-mechanical
properties remains unknown to date. Based on surveys of other CCKBOs (Lellouch
et al., 2013) only bulk thermal inertias probing diurnal skin depths are
reasonably known. For the solutions displayed in Section we have adopted a
representative value of ${\cal I}\approx 2.5$ tiu, but we keep in mind that
these values may be larger based on the discussion found in section 6.1.

Nevertheless, these thermal solutions with relatively low values of ${\cal I}$
might suggest something about Arrokoth’s erosive character. Erosion driven by
thermal cycling – wherein cracks nucleate and grow due to a material’s non-
zero coefficient of thermal expansion – might be a relevant effect on Arrokoth
over the course of its natural history. Such a process has been hypothesized
in order to explain the nature and origin of erosion on for cometary bodies
(e.g., for 67P/Churyumov-Gerasimenko, Attree et al., 2018).

We focus here on thermomechanical erosion and do not consider thermal erosion
driven by sublimation since volatile ices are neither known to be present on
Arrokoth’s surface nor would be expected to be extant there (see previous
section). Among other things, however, estimating actual rates of erosion
requires knowing what is the local temperature temporal profile, the materials
coefficient of thermal expansion, its effective porosity, Young’s modulus,
Poisson’s ratio, and if the thermal cycling is on a long enough timescale, its
viscous relaxation properties like the materials activation energy (Mellon,
1997; Molaro et al., 2015, see latter reference for a recent comprehensive
discussion). Nevertheless, we speculate upon the relative amount of thermally
driven mechanical erosion _over the course of one orbital time_ assuming the
near surface of the body is uniformly composed. Fig. 11 indicates that the
highest minimum temperature ($\sim 20-30$K) over the course of one orbit
occurs along a narrow line of surface located within the tropical zone and,
moreover, appears to trace along the rims of putative craters or depressions
found in that region. The implied extreme temperature variation is shown in
Fig. 13, which shows for each facet the difference of the maximum and minimum
temperatures achieved over the course of one orbit, i.e.,

$\Delta T_{{}_{i}}\equiv{\rm max}(T_{{}_{i}})-{\rm min}(T_{{}_{i}}).$ (42)

The figure indicates that the bulk of the polar regions experiences $\Delta
T_{{}_{i}}$ close to 50 K, while the tropical zones experience $\Delta
T_{{}_{i}}$ closer to 30K. The aforementioned narrow line within the tropical
zone has a temperature variation of about only 10K. With all other variables
held fixed, it stands to reason that the thermal erosion is greatest across
the polar regions where $\Delta T_{{}_{i}}$ is large, and significantly
weakens upon entry into the tropical zone and even further weakens as one
approaches the narrow line of real estate within in which $\Delta
T_{{}_{i}}\approx 15-20$K. While this relative prediction is reasonable, the
actual amount of erosion that has taken place over Arrokoth’s natural history
cannot be ascertained until the actual composition of its surface materials
are unambiguously determined, together with information about their thermal
and viscoelastic properties, with particular attention to its rheology.

As a final remark, we note that the thermal driving is likely most extreme on
Arrokoth’s equatorial latitudes during its equinoxes since these regions are
getting extreme illumination variations on the the short 15.9 hour diurnal
period $P_{{}_{{\rm day}}}$ ($\omega_{{}_{{\rm day}}}=2\pi/P_{{}_{{\rm
day}}}$). But the thermomechanically driven erosion taking place on this short
timescale is also restricted to the upper 1-2mm of surface since the thermal
wave lengthscale $\ell_{{}_{{\rm day}}}$ is shortened by comparison to
$\ell_{{}_{{\rm orb}}}$ ($\approx 1$m) by a factor of $\sqrt{\omega_{{}_{{\rm
orb}}}/\omega_{{}_{{\rm day}}}}\approx 1/400$.

###  6.5 Other Various Caveats

Despite the intrinsic scientific value of these solutions, several aspects of
the thermal model method utilized here could be improvement. We acknowledge a
few here with the aims of developing even more refined solutions in the
future. We have assumed the thermal inertia (and by implication, the thermal
conductivity) is a constant independent of temperature and porosity. Indeed
the time-asymptotic method we have employed relies on this feature, which
permits using fast Fourier transform techniques to develop solutions (Titus &
Cushing, 2012; Schloerb et al., 2015; White et al., 2016). However, although
not done here, a temperature dependence of the thermal inertia may be formally
built into the method at the cost of added iterations to the solution process
– a matter we will develop further with future follow-up work. We have also
assumed constant stratigraphy and/or subsurface porosity, which could also be
unrealistic. But similar to our point above this can be similarly remedied by
iteration or by directly building stratigraphy into the solution itself. We
demonstrate how this is done in Appendix C where we develop a subsurface
thermal solution for a material composed of two different conductivities where
the transition occurs $h$ meters below the surface. The solution may be
developed analytically, as was done for the discussion in the text, and coded
directly into the thermal solver. This procedure naturally generalizes to
multiple layers with differing properties.

####  6.5.1 On the use of a single albedo value

New Horizons imaging showed the existence of several bright rings and patches
painting its surface (e.g., see the bright patches in Fig. 4, also see Spencer
et al., 2020), with normal reflectances of these bright patches being about
twice that of the rest of the terrain. Given that these images were taken at
high solar phase angle, the observed relatively enhanced brightness in these
regions is likely “the result of the extreme illumination conditions rather
than the intrinsic albedo of the surface,” (Stern et al., 2021). In light of
this and in addition to determining Arrokoth’s mean hemispherical albedo,
Hofgartner et al. (2021) have produced an albedo map of the closest approach
image (see their Fig. 4) in which the bright patches are assessed to have
albedos as high as 0.09. Provided all other thermophysical properties are the
same across bright and dark patches, we find that the temperatures in these
bright regions are overestimated by less than 0.5K: if $T_{{\rm dark}}$ is the
estimated temperature of the bright terrain assuming the lower albedo
$(A_{{\rm dark}}=0.063$) and $T_{{\rm bright}}\equiv T_{{\rm dark}}-\delta T$
would be its true correct value based on $(A_{{\rm bright}}=0.09)$, where
$\delta T$ is the correction factor, then

$\left(\frac{T_{{\rm dark}}-\delta T}{T_{{\rm
dark}}}\right)^{4}=\frac{1-A_{{\rm dark}}}{1-A_{{\rm bright}}},$ (43)

predicting $\delta T\big{/}T_{{\rm dark}}\approx 0.007$, which for $T_{{\rm
dark}}\approx 60$K amounts to $\delta T=0.42$K.

####  6.5.2 On the neglect of reflected light

In our calculations we have included thermal IR in the re-radiation
calculation and neglected reflected light. This neglect has minor consequences
on the final determined temperatures following the same line of thinking as
above, where we find the error is even weaker. If we treat the reflected light
as a Lambertian process, then the light received at a given facet will be
traced by the IR re-radiated light from all other facets. We showed earlier
that for a given facet $i$ the energy received from all other facets
($f_{i,{\rm rr}}$) amounts to about 10 percent of the locally absorbed solar
insolation ($f_{{}_{i,\odot}}$), i.e., $f_{i,{\rm rr}}\approx\chi
f_{{}_{i,\odot}}$, where $\chi\leq 0.1$. We can estimated the amount of
reflected light neglected in this accounting to be about $Af_{i,{\rm rr}}$,
such that the correct amount of light received from the surrounding landscape
is better estimated by $\tilde{f}_{i,{\rm rr}}\approx f_{i,{\rm
rr}}\big{/}(1-A)$. Assuming that all other thermophysical properties remain
the same, the correct temperature $T_{i,{\rm cor}}=T_{i}+\delta T$ compares to
the estimated temperature $T_{i}$ via the same Stefan-Boltzmann law argument
applied in the previous section with

$\left(\frac{T_{i}+\delta T_{i}}{T_{i}}\right)^{4}=\frac{f_{i,\odot}+f_{i,{\rm
rr}}}{f_{i,\odot}+\tilde{f}_{i,{\rm
rr}}}=\displaystyle\frac{1+\chi}{\displaystyle{1+\frac{\chi}{1-A}}},$ (44)

which, given our estimate for $\chi$, means $\delta T_{i}/T_{i}\approx
0.0016$, which for $T_{i}\approx 60K$, implies that neglecting reflected light
on Arrokoth amounts to underestimating the local temperature by $\delta
T_{i}\approx 0.1$K. In summary it is conceivable that Arrokoth is close to a
blackbody, i.e., one that is marginally reflective and radiating most of its
heat in IR.

####  6.5.3 On the expiration of surface N2 ice on Arrokoth

We can ask about the fate of volatile ices like N2 and CO (with similar
properties) on an airless body like Arrokoth following some of the thinking
used to address similar concerns about the nature of Oumuamua (Jackson &
Desch, 2021). With typical solar insolation values based on Arrokoth’s
semimajor axis, i.e., $\tilde{f}_{\odot}\approx$ 0.7 W/m2, a vacuum exposed
surface layer of N2 ice on Arrokoth would sublimate at a rate of
3.5m/m2/orbit. This conservative figure derives from assuming an N2 ice albedo
of $A=$0.9 (i.e., based on Pluto observations, Grundy et al., 2016; Schmitt et
al., 2017) and an enthalpy of N2 sublimation ${\cal L}=225$kJ/kg that, if we
assume a surface layer of N2 ice converts all of its received insolation into
sublimation (and adjusting its temperature accordingly), would result in a
mass flux $\dot{\Sigma}=(1-A)\tilde{f}_{\odot}/{\cal L}=3.4\times
10^{-7}$kg/m2/s $\approx 3.24\times 10^{3}$kg/m2/Arrokoth orbit. Under these
vacuum conditions a sublimating block of ice must satisfy the relationship
between mass-flux into a vacuum and a material’s vapor pressure,$P_{{\rm
vap}}(T)$, i.e.,

$\dot{\Sigma}v_{k}=P_{{\rm vap}}(T),$ (45)

(e.g., Lebofsky, 1975), where $v_{k}$ is the typical Boltzmann velocity. An
estimate based on compiled vapor pressure data for N2 (Fray & Schmitt, 2009)
shows that this balanced state is satisfied for $T\approx 27$K. Based on the
shape model utilized here Arrokoth’s total surface $A_{{\rm surf}}\approx
1406$km2. Thus, if Arrokoth were uniformly covered with N2 ice, then it would
be losing a total of about 481 kg of N2 per second $\leftrightarrow 1.6\times
10^{28}$ particles per second, well above the mass loss rate upper limits
based on New Horizons’ observation of the body ($<10^{24}$particles/s, Stern
et al., 2019). These considerations preclude the plausibility that N2 and CO
can remain on the surface of Arrokoth for any appreciable amount of time. CH4,
which is relatively less volatile by a factor of $10^{3}$ (based on its vapor
pressure behavior at similar temperatures) would also result in particle flux
rates exceeding those observed by New Horizons. This analysis says nothing
about the plausibility for reservoirs of volatile ices deep within Arrokoth.
Addressing this matter requires further analysis.

##  7 Summary

We have developed a temperature model of Arrokoth based on the $10^{5}$ facet
model published in Spencer et al. (2020). For the bulk of the solutions
displayed here, especially in Section 4, we have assumed a single value of
${\cal I}\approx 2.5$tiu, together with a infrared emissivity of $0.9$ and an
albedo of 0.06. The choice of ${\cal I}\approx 2.5$ is consistent with the low
end of thermal inertia values for amorphouse H2O ice as detailed at length in
section 6.1. Although direct temperature measurements were not possible with
New Horizons, for the day of the closest encounter we predict that the day
side approach hemisphere surface temperatures were in the range of 57-60K.
According to the simple model of section 2 (e.g., Fig. 2) the predicted
daytime (“summer time”) surface temperature has a very weak dependence in the
thermal inertia range of 1 tiu $<{\cal I}<$ 10 tiu. We thus expect this quoted
daytime temperature range to hold for the full range of likely ${\cal I}$
values. On the other hand, for the specific ${\cal I}=2.5$ tiu value we
predict that the obscured winter side of Arrokoth had surface temperatures in
the range of 10-15K, but following the same above reasoning, based on the
results shown in Fig. 2 together with the likely range for ${\cal I}$, we
predict that the night side surface temperature could lie within the slightly
wider range of values of 10-20K.

Arrokoth’s tropical zones receive the lowest yearly averaged insolation but
exhibits the highest orbitally averaged temperatures compared to Arrokoth’s
more extreme polar regions. In fact, the extreme temperature variations
experienced by facets in Arrokoth’s polar zones ($\Delta T\approx 50$K)
greatly exceed the corresponding variations in its equatorial zones ($\Delta
T\approx 10-30$K). The property of our thermal model leads us to conjecture
that the amount of thermal erosion in the tropical zones should be
significantly lower than the erosion suffered in the polar regions.

We also find that the amount of energy received on a typical facet that is
sourced from surface reradiation typically amounts to about 5% of the total
irradiation, although it should be kept in mind – as a cross comparative
inspection of the histograms of Fig. 7 shows – there is a spread in this
figure from facet to facet, which ultimately depends upon how much of a given
facet is obscured by the rest of the body. A small adopted value of ${\cal I}$
means that Arrokoth’s Spencer Number is also small ($\approx 0.01$). In other
words, this highly insulating body has very little downward directed thermal
conduction from radiation received during total direct illumination (and
subsequently returning it at night), amounting to about 0.5% of the typically
received insolation flux budget.

Based on a generalized application of our thermal modeling we find that the
brightness temperature $T_{{\rm b,obs}}=29\pm 5$K measured by New Horizons’
REX instrument (Grundy et al., 2020, also see Bird et al., 2022) is consistent
with porous $(0.6<p<0.8$) subsurface amorphous H2O ice – possibly coated with
tholins – within the thermal inertia range $1\ {\rm tiu}<{\cal I}<$10-20 tiu
for given X-band emissivity $\big{<}E_{{\rm eff}}\big{>}=1.0$. For X-band
emissivity $\big{<}E_{{\rm eff}}\big{>}=0.9$, which is similar to the low
albedo (and possibly captured) KBO Phoebe, our thermal modeling predicts a
wider range of permissible thermal inertias, with $1\ {\rm tiu}<{\cal I}<$ 300
tiu. However the high-end of the ${\cal I}$ ($>50$ tiu) values are disfavored
for two reasons: (i) it would imply that either very large (1-10m) or very
small (sub-micron) amorphous H2O ice grains are dominant constituents of
Arrokoth’s near surface – sizes which are hard to reconcile from an origins
planetesimal formation perspective – and, (ii) it might imply that the near
surface H2O ice grains are in cubic crystalline form, but this would mean that
Arrokoth’s surface experienced sustained periods of relatively elevated
temperatures exceeding 70K (corresponding to the amorphous-cubic transition
occurring on times shorter than Arrokoth’s natural age). Such high thermal
inertias might be characteristic of low-porosity tholins, but that is
currently unknown and remains to be more completely studied under relevant
laboratory conditions.

As such we favor the interpretation that the REX $T_{{\rm b,obs}}$ measurement
indicates that Arrokoth harbours porous hydrocarbon-coated amorphous H2O ice
grains of rough size $\sim 0.1-1$mm, with $1\ {\rm tiu}<{\cal I}<$10-20 tiu,
and characterized by an X-band emissivity in the range 0.9 and 1.

Acknowledgements. O.M.U. and the rest of the New Horizons team acknowledge
support from the New Horizons Kuiper Belt Extended Mission grant for support
in producing this work. J.T.K. acknowledges that portion of the research was
carried out at the Jet Propulsion Laboratory, California Institute of
Technology, under a contract with the National Aeronautics and Space
Administration (80NM0018D0004). J.T.K. also acknowledges supported by the
California Institute of Technology Joint Center for Planetary Astronomy
postdoctoral fellowship. O.M.U. acknowledges fruitful conversations with A.
Jindal and A. Hayes (Cornell University). We are indebted to insights shared
with us by B. Butler (NRAO) and M.A. Gurwell (CFA). The authors acknowledge
and are grateful for the helpful comments made by the two reviewers of the
manuscript.

##  Appendix A Details of Shape Analysis of Section 3.2

This section is concerned with determining how much a given piece of landscape
is visible to the rest of the landscape (a “who-sees-who” list), and is meant
to address how the integral expression in Eq. (15) is calculated. The
considerations here reference Fig. 22. Each facet $i$ has an outwardly
pointing normal $\hat{\bf n}_{{}_{i}}$ and associated surface area
$s_{{}_{i}}$ such that $d\hat{\bf s}_{i}=s_{{}_{i}}\hat{\bf n}_{{}_{i}}$. Our
goal is to determine which facets $j$ are visible to facet $i$. This involves
ruling out all (1) facets $j$ whose unit normals $\hat{\bf n}_{{}_{j}}$ face
away from facet $i$, (2) rejecting all facets $j$ that are below the facet
$i$’s horizon, and finally (3) discarding all facets $j$ that are blocked by
intervening facets. We treat each of these in order:

We first determine a separation vector ${\bf r}_{{}_{ij}}={\bf
r}_{{}_{j}}-{\bf r}_{{}_{i}}$, and further parse this in terms of a scalar and
unit vector ${\bf r}_{{}_{ij}}=r_{{}_{ij}}\hat{\bf m}_{{}_{ij}}$, in which
$r_{{}_{ij}}=|{\bf r}_{{}_{ij}}|$ and $\hat{\bf m}_{{}_{ij}}$ is a unit vector
pointing from the center of facet $i$ to the center of facet $j$: (1) a facet
$j$ is deemed below the horizon if $\cos\theta_{{}_{i}}=\hat{\bf
m}_{{}_{ij}}\cdot\hat{\bf n}_{{}_{i}}<0$, (2) a facet $j$ is pointed away if
$\hat{\bf m}_{{}_{ij}}\cdot\hat{\bf n}_{{}_{j}}>0$, for which we define
$\cos\theta_{{}_{j}}\equiv-\hat{\bf m}_{{}_{ij}}\cdot\hat{\bf n}_{{}_{j}}$.
Therefore in order for a facet $j$ to be a candidate for visibility by facet
$i$, both $\cos\theta_{{}_{i}}>0$ and $\cos\theta_{{}_{j}}>0$. It then follows
that all candidate facets $j$ passing the above two criteria are tested for
actual visibility using a ray-tracing routine. The prototype code we have
built is run on Matlab 2014a and we have therefore used a community wide
tested ray-tracing Matlab routine called inpolyhedron, which ascertains
whether not elements of a given ray of points are inside a closed surface.

For all the facets $j$ visible to facet $i$, we determine the local altitude
of the highest visible facet by identifying for which facet $j$ is
$\cos\theta_{i}$ maximum, and we assign the value of this special altitude
$\beta_{{}_{i}}=\pi/2-\theta_{i}({\rm max})$. This information is used in
speeding up the solar insolation calculation of the next section. In some
notable instances a given facet $i$ sees no other facets of the body, in which
case its $\beta_{{}_{i}}=0$.

The luminosity apparent to facet $i$ of a _visible_ flat Lambertian facet $j$,
with area $s_{{}_{j}}$ and inclined at an angle $\cos\theta_{j}$, emitting
with a flux $f_{{}_{j}}$ is $f_{{}_{j}}s_{{}_{j}}\cos\theta_{j}$. Defining
$F_{{}_{ij}}$ as the flux absorbed by facet $i$ from facet $j$, it follows
that

$F_{{}_{ij}}=f_{{}_{j}}{\rm K}_{{}_{ij}},\qquad{\rm
K}_{{}_{ij}}=\frac{s_{{}_{j}}\cos\theta_{j}\cos\theta_{i}}{2\pi
r_{{}_{ij}}^{2}},$ (A1)

where the factor $2\pi$ accounts for the fact that emitted radiation spreads
out only across 2$\pi$ steradians. We note that the expression in Eq. (A1)
likely badly quantifies the received irradiation from (say) neighboring facets
especially if the angles formed between their planes of intersection are
highly acute. This problem would be especially acute if the quality of the
local topography _on the scales of interest_ appeared to be fractal. This is
not the case here. By identifying $f_{{}_{j}}=\varepsilon\sigma
T_{{}_{j}}^{4}$, then it also follows that the body integral term in Eq. (15)
becomes

$\int_{\partial S}f_{{}_{j}}S_{{}_{ij}}d\hat{\bf
s}_{j}\rightarrow\sum_{\forall j\ {\rm visible}}f_{{}_{j}}{\rm K}_{{}_{ij}},$
(A2)

where, by implication, $S_{{}_{ij}}d\hat{\bf s}_{j}\longrightarrow{\rm
K}_{{}_{ij}}$. In practice ${\rm K}_{{}_{ij}}$ is constructed as a sparse
matrix in which the column elements $j$ are non-zero for those facets that are
visible to $i$. By its construction a given facet is not visible to itself,
therefore ${\rm K}_{{}_{ii}}=0$. Finally, we note that although this analysis
is expensive ($\sim 5$ hours when run as a parallel computation on a 4-core
3.1 GHz Intel system), it is only done once for a given _static_ shape model.
All other calculations, including parameter studies and so forth that are done
hereafter, utilize ${\rm K}_{{}_{ij}}$.

In Fig. 3 we show the results of the shape analysis based on the 105 faceted
shape model described in the previous section. The figure depicts a facet-by-
facet map of the percent of open sky self-obscured by Arrokoth. The percent is
based on the fraction of $2\pi$ steradians occupied by other surface facets.
We see how the routine identifies craters as places of relatively high
obscuration ($\sim\%30$). The neck region on the well-imaged approach
trajectory (top panel) also shows similar degrees of self-obscuration.
However, the neck region on the departure trajectory, which was not directly
imaged and whose structure is inferred from model fitting procedures described
in Spencer et al. (2020) and is probably less reliable as a result, is
identified as having much greater obscuration in the vicinity of $\sim\%60$ or
so. It is not a surprise given Arrokoth’s flattened structure that the
relative amount of self-obscuration is generally very weak (also see the
distribution of these percentages in the first panel of Fig. 12).

Figure 22: This diagram describes the shape analysis discussed in Section 3.2.
The shape model depicted here is the low order 2000 facet model analyzed in
Grundy et al. (2020). However, the methods described herein are applied to the
$10^{5}$ facet model in Spencer et al. (2020).

##  Appendix B An example calculation using the Fourier method

We detail the Fourier method solution procedure on a relatively simplified
problem. We suppose we have a 1 dimensional layer with the same setup as
described in the text. We imagine that the surface ($z=0$) is illuminated by a
function $f$ with period $2\pi/\omega$. The surface material’s emissivity is
1. The interior temperature satisfies the heat equation

$\rho C_{p}\partial_{t}T=K\partial_{z}^{2}T,$ (B1)

The boundary conditions at the top and bottom (respectively) are

$\sigma
T^{4}\Big{|}_{z=0}+K\partial_{{}_{z}}T\big{|}_{z=0}=f,\qquad\partial_{z}T\big{|}_{z\rightarrow\infty}=0.$
(B2)

We assume all other coefficients in the above equations are constant. We
represent the interior solution by a Fourier series

$T=T_{{}_{0}}+\sum_{n=1}^{\infty}\Theta_{{}_{n}}e^{k_{n}z-i\omega nt}+{\rm
c.c.}.,\qquad k_{{}_{n}}=\sqrt{\rho C_{p}i\omega n\Big{/}K}.$ (B3)

The above solution satisfies the heat equation for all points $z\leq 0$,
including the bottom boundary condition. What remains to be determined are the
coefficients $\Theta_{{}_{n}}$, which can be solved for by satisfying the
surface boundary condition. We approach this in the following way. The
vertical thermal gradient at $z=0$ is given by

$K\partial_{{}_{z}}T\big{|}_{z=0}=K\sum_{n=1}^{\infty}k_{{}_{n}}\Theta_{{}_{n}}e^{-i\omega
nt}++{\rm c.c.}.$ (B4)

We can decompose the insolation forcing also as an infinite sum of Fourier
harmonics, that is to say

$f=f_{{}_{0}}+\sum_{n=1}^{\infty}f_{{}_{n}}e^{-i\omega nt}+{\rm c.c.}.$ (B5)

The top boundary condition therefore becomes

$\sigma
T^{4}\Big{|}_{z=0}-f_{0}+\sum_{n=1}^{\infty}(Kk_{{}_{n}}\Theta_{{}_{n}}-f_{{}_{n}})e^{-i\omega
nt}+{\rm c.c.}.=0,$ (B6)

A complete solution is had once the countably infinite number of complex
valued coefficients $\Theta_{{}_{n}}$ are determined. In practice one assumes
a maximum number $N$ of these coefficients for a sufficiently converged
solution, which is assessed as such once successive changes in the final state
meets some requisite set of minimum error criteria. Thus, the solution is
determined for the set of $\Theta_{{}_{n}}$ that solve Eq. (B6). Identifying
the function $F(T)\equiv\sigma T^{4}\Big{|}_{z=0}$, together with its Fourier
decomposition as

$F=F_{{}_{0}}+\sum_{n=1}^{\infty}F_{{}_{n}}e^{-i\omega nt}+{\rm c.c.}.$ (B7)

we find that Eq. (B6) becomes after truncating up to $N$ sinusoids

$F_{{}_{0}}-f_{0}+\sum_{n=1}^{N}(F_{{}_{n}}+Kk_{{}_{n}}\Theta_{{}_{n}}-f_{{}_{n}})e^{-i\omega
nt}+{\rm c.c.}.=0.$ (B8)

Taking the (fast) Fourier transform of the above yields a set of $N$ nonlinear
algebraic equations that must be simultaneously solved, namely

$F_{{}_{n}}(\Theta_{{}_{m}})+Kk_{{}_{n}}\Theta_{{}_{n}}-f_{{}_{n}}=0,\qquad
k_{{}_{0}}=0,$ (B9)

for each value of $n=0,1,\cdots N$. The term
$F_{{}_{n}}=F_{{}_{n}}(\Theta_{{}_{m}})$ is a nonlinear function of all the
unknown coefficients $\Theta_{{}_{m}}$. As such, the above set of N-coupled
nonlinear equations must be solved using a Newton-Raphson (or gradient
descent) type of nonlinear equation solver. In practice, this is done
efficiently by assuming an initial guess for $\Theta_{{}_{n}}$, then backward
transforming $T(z=0)$ from frequency space into temporal space, followed by
calculating $F(T)=\sigma T^{4}\big{|}_{z=0}$, then taking $F$ and forward
(fast) Fourier transforming it into frequency space to determine each
$F_{{}_{n}}$, followed by assessing a solution to Eq. (B9). This usually
involves an iterative procedure which requires re-executing the same set of
forward-backward Fourier transform procedures to determine an improved
solution, until error tolerances are met. This solution procedure is easily
generalizable to nonlinear values of the conductivity coefficient $K$, say,
which are dependent on the values of $T$.

##  Appendix C Layered thermal solution

We suppose that the surface is composed of two layers with a transition
occurring at $z=-h$ (where $h>0$). For the sake of this illustration let us
suppose that the only difference between the layers is the thermal
conductivity. In keeping with our convention in the text, we set the
conductivity to be $K$ for $-h<z<0$, and $K^{(-)}$ for $z<-h$. We therefore
divide the solution to the heat equation Eq. (13) into ones appropriate to
these two levels:

$\Theta_{i}=\left\\{\begin{array}[]{cc}\sum_{n=0}^{N}\left(\Theta_{i,n}e^{k_{n}z}+\Xi_{i,n}e^{-k_{n}z}\right)e^{-i\omega
nt}+{\rm c.c.},&-h<z\leq 0\\\
\sum_{n=0}^{N}\Theta_{{}_{i,n}}^{(-)}e^{k_{n}^{(-)}(z+h)}e^{-i\omega nt}+{\rm
c.c.},&z\leq-h\end{array}\right.,$ (C1)

where $\Xi_{i,n}$ and $\Theta_{{}_{i,n}}^{(-)}$ are new frequency amplitudes
that must be determined in terms of $\Theta_{i,n},h$ and $K^{(-)}$. The lower
layer wavenumber $k_{n}^{(-)}$ is analogous to $k_{n}$, as defined in Eq.
(23), where

$k_{n}^{(-)}=\left(\frac{i\rho C_{p}\omega n}{K^{(-)}}\right)^{1/2},\qquad{\rm
with}\ \ {\rm Re}\left(k_{n}^{(-)}\right)>0.$ (C2)

The thermal solution form of Eq. (C1) is entirely analogous to Eq. (22). The
solution for the lowest layer is designed to have zero flux as
$z\rightarrow-\infty$. The solutions across the layer separation $z=-h$ are
connected to one another for each value of $n$ by imposing the continuity of
temperature across the level, i.e.,

$\Theta_{i,n}e^{-k_{n}h}+\Xi_{i,n}e^{k_{n}h}=\Theta_{{}_{i,n}}^{(-)},$ (C3)

as well as the continuity of thermal flux across the layer, i.e.

$\Big{[}K\partial_{z}T\Big{]}^{z\rightarrow-h^{+}}_{z\rightarrow-h^{-}}=0,$
(C4)

which means that for each index $n$ we require

$k_{n}\Theta_{i,n}e^{-k_{n}h}-k_{n}\Xi_{i,n}e^{k_{n}h}=k_{n}^{(-)}\Theta_{{}_{i,n}}^{(-)}.$
(C5)

Solving these two simultaneous equations admits

$\Xi_{i,n}=\Theta_{i,n}\frac{1-\sqrt{K/K^{(-)}}}{1+\sqrt{K/K^{(-)}}}e^{-2hk_{n}},\qquad\Theta_{{}_{i,n}}^{(-)}=\displaystyle{\frac{2\Theta_{i,n}e^{-hk_{n}}}{1+\sqrt{K/K^{(-)}}}}.$
(C6)

Therefore, the temperature at the surface is

$T_{{}_{i}}=\sum_{n=0}^{N}\Theta_{i,n}\left(1+\frac{1-\sqrt{K/K^{(-)}}}{1+\sqrt{K/K^{(-)}}}e^{-2hk_{n}}\right)e^{-i\omega
nt}+{\rm c.c.},$ (C7)

while the conductive flux is

$\Lambda_{{}_{i}}=K\partial_{z}T_{{}_{i}}=\sum_{n=1}^{N}k_{n}\Theta_{i,n}\left(1-\frac{1-\sqrt{K/K^{(-)}}}{1+\sqrt{K/K^{(-)}}}e^{-2hk_{n}}\right)e^{-i\omega
nt}+{\rm c.c.}$ (C8)

Consequently, the solution method described in Sec. 3.4 may be executed
exactly as it is there described, except for replacing Eq. (24) with Eq. (C8),
and the surface temperature $T_{{}_{i}}$ replaced by the expression in Eq.
(C7). In reference to the illustrative case described in Appendix B, Eq. (B9)
would be rewritten instead as

$F_{{}_{n}}(\Theta_{{}_{m}})+Kk_{n}\Theta_{n}\left(1-\frac{1-\sqrt{K/K^{(-)}}}{1+\sqrt{K/K^{(-)}}}e^{-2hk_{n}}\right)-f_{{}_{n}}=0,\qquad
k_{{}_{0}}=0.$ (C9)

The solution method described here may be straightforwardly generalized to
multiple layers. The coefficients may be calculated using any symbolic
mathematics software or using any prepackaged linear system solver.

##  Appendix D Subsurface Radiative Transfer Details

In order to connect these interior thermal solutions to the observed
brightness temperature, $T_{b,{\rm obs}}=29\pm 5$K, measured by New Horizons
REX instrument, we have implemented the radiative transfer calculation
detailed at length in de Kleer et al. (2021). A brightness temperature $T_{b}$
is defined as the temperature of a blackbody that produces the observed
radiative intensity at frequency $\nu$. We designate the Arrokoth disk
averaged REX observed intensity by ${I}({\rm REX})$, which is the sum of the
disk-averaged Arrokoth emitted intensity at frequency $\nu$
($\tilde{I}_{\nu}$) plus the intensity received from the cosmic microwave
background ($I_{\nu,{\rm cmb}}$). The identification between ${I}({\rm REX})$
and $T_{b,{\rm obs}}$ comes from the solution of

$\displaystyle{\cal I}({\rm REX})=\tilde{I}_{\nu}+I_{\nu,{\rm cmb}}$
$\displaystyle=$ $\displaystyle\frac{\displaystyle
2h\nu^{3}\big{/}c^{2}}{\displaystyle{\exp\left[{{\frac{h\nu}{kT_{b,{\rm
obs}}}}}\right]-1}}+\frac{\displaystyle
2h\nu^{3}\big{/}c^{2}}{\displaystyle{\exp\left[{{\frac{h\nu}{kT_{{\rm
cmb}}}}}\right]-1}},$ (D1)

where $h$ is Planck’s constant, $c$ speed of light and where the frequency
$\nu=c/\lambda$. Given that the X-band REX observation occurs at a wavelength
$\lambda=0.04$m, places this emission in the Rayleigh-Jeans end of the
blackbody spectrum as $T_{\nu}\equiv h\nu/k\approx 0.36$K $\ll T_{b,{\rm
obs}}$ and $T_{{\rm cmb}}$, which simplifies the above into

${I}({\rm REX})=\frac{2k\nu^{2}}{c^{2}}\Big{(}T_{b,{\rm obs}}+T_{{\rm
cmb}}\Big{)}.$ (D2)

$\tilde{I}_{\nu}$ is the average of the emitted intensity
($\tilde{I}_{\nu,i}$) across the set of all spacecraft-viewable facets $i$
_weighted_ by each facet’s spacecraft projected area $\sigma_{i}$. We do not
include any Gaussian weighting that would take into account of the REX beam.
For each facet $i$ with unit normal vector ${\bf\hat{n}}_{i}$ we define
$\theta_{i}$ as the emission angle of the rays emitted from the surface that
reach the spacecraft, from which it follows that if each facet has total area
$s_{i}$, then $\sigma_{i}=s_{i}\cos\theta_{i}$. For geometrical reference see
both Fig. 23 as well as Fig. 22 in Appendix A. By construction it follows that

$\tilde{I}_{\nu}=\frac{\displaystyle\sum_{\forall i\ {\rm
visible}}\tilde{I}_{\nu,i}\sigma_{i}}{\displaystyle\sum_{\forall i\ {\rm
visible}}\sigma_{i}},$ (D3)

where all facets viewable by the spacecraft at the time of observations is
derived directly from the network analysis that went into developing
Arrokoth’s thermal solutions as described in Section 3 and Appendix A.

Figure 23: Diagram for simple radiative transfer calculation.