# Constraining Thermal Emission of Pluto’s Haze From Infrared Rotational
Lightcurves

Linfeng Wan (万霖丰)  Department of Earth and Planetary Sciences, University of
California Santa Cruz, CA 95064, USA  Xi Zhang  Department of Earth and
Planetary Sciences, University of California Santa Cruz, CA 95064, USA  Jason
D. Hofgartner  Southwest Research Institute, Boulder, CO 80302, USA

###### Abstract

The rotational lightcurves of the Pluto-Charon system were previously believed
to be solely attributed to their surfaces. However, a proposed scenario of
haze cooling (Zhang et al., 2017) suggests that the atmospheric haze of Pluto
could significantly contribute to mid-infrared emission, which calls for a
revisit of previous analyses. In this study, we employ a Bayesian retrieval
approach to constrain the haze emission from the rotational lightcurves of the
Pluto-Charon system. The lightcurves were observed by the Spitzer and Herschel
telescopes at 24 and 70 $\mu$m, and were combined with the latest surface
albedo maps of Pluto and Charon from the New Horizons spacecraft. Our results
show that including the haze emission is consistent with all current
observations, with the best-fit haze flux around 1.63 mJy. This is in
agreement with the composition of Titan-like tholins. However, the “surface
only” scenario, which excludes the haze contribution, can still explain the
observations. We conclude that the current data at 24 $\mu$m cannot constrain
Pluto’s haze emission due to the degeneracy with Charon’s surface emission.
Regardless, some surface properties of Pluto are well constrained by the shape
of the lightcurves, with a thermal inertia of approximately 8–10 MKS and a
relatively low CH4 emissivity of 0.3–0.5. We suggest that observations by the
JWST telescope at 18 $\mu$m, which can resolve Pluto from Charon, could
directly probe the haze emission of Pluto due to the low surface emission at
that wavelength.

planets: surface – planets: Pluto – techniques: retrieval

##  1 Introduction

Rotational lightcurves provide insights into the surface properties and
inhomogeneities of Pluto and Charon as the surface units rotate in and out of
view. While reflection lightcurves mainly reveal the body’s albedo and
scattering properties, emission lightcurves at infrared and radio wavelengths
offer valuable information about the surface and subsurface, such as
temperature, thermal inertia, and emissivity. Before the New Horizons flyby,
Pluto’s surface albedo distribution could only be resolved with poor precision
using techniques such as Pluto-Charon mutual events (e.g., Buie et al., 1992;
Young & Binzel, 1993; Young et al., 1999, 2001) and direct imaging by HST
(e.g., Buie et al., 1997; Stern et al., 1997; Buie et al., 2010a, b).

Nevertheless, by combining reflection and emission lightcurves from HST, ISO,
Spitzer, and Herschel, Lellouch et al. (2000, 2011, 2016) found that Pluto’s
surface temperature is non-uniform. They also constrained the thermal inertia
and emissivity of the surfaces of Pluto and Charon. Multiple solutions have
been identified, primarily due to the inability to resolve the Pluto-Charon
system and the various choices of Pluto’s surface map (Figure 4 in Lellouch et
al. 2011). Although radio telescopes such as the Atacama Large
Millimeter/submillimeter Array (ALMA), Submillimeter Array (SMA), and Karl G.
Jansky Very Large Array (VLA) are able to observe Pluto and Charon
individually, the seasonal effect and the mixing of subsurface and surface
information at radio wavelengths have complicated data interpretation
(Lellouch et al., 2016).

The New Horizons flyby of the Pluto-Charon system in 2015 provided a wealth of
new information that calls for revisiting previous rotational lightcurve data.
High-definition images of Pluto and Charon obtained by the Long Range
Reconnaissance Imager (LORRI, Cheng et al., 2008) and the Multispectral
Visible Imaging Camera (MVIC, Reuter et al., 2008) provided detailed spatial
distribution of geological units, and further spectral analysis of data by the
Linear Etalon Imaging Spectral Array (LEISA, Reuter et al., 2008) revealed
complex ice compositions on different surface units composed of N2, CO, CH4,
H2O, etc. It was found that Pluto’s surface has complex volatile distributions
that are closely related to ice sublimation, gas condensation, and glacial
flow (Grundy et al., 2016; Moore et al., 2016). On the other hand, Charon has
a relatively uniform surface with a reddish northern polar region (Grundy et
al., 2016; Stern et al., 2017). The preliminary surface albedo maps of Pluto
and Charon based on New Horizons data, derived by Buratti et al. (2017, 2019),
contain much greater detail than the surface maps in previous studies.
Therefore, it is important to revisit previous thermal rotational lightcurve
data based on the new albedo maps, which is the first motivation of this
study.

The second motivation is related to the hazy atmosphere of Pluto.
Traditionally, the rotational lightcurve analysis of the Pluto-Charon system
was focused only on their surfaces. However, during the New Horizons flyby, it
was discovered that Pluto’s atmospheric temperature was much colder than
previously thought based on the gas-only model. Proposed solutions include
cooling from the high concentration of super-saturated water vapor (Strobel &
Zhu, 2017) or from the atmospheric haze (Zhang et al., 2017; Wan et al., 2021)
that was clearly visible in the New Horizons images (Gladstone et al., 2016).
If haze is the main coolant in Pluto’s atmosphere, it is predicted that the
thermal emission from the haze could significantly contribute to the observed
mid-infrared flux, such as that observed at 23.68 (hereafter 24) $\mu$m by
Spitzer (Lellouch et al., 2011), while the haze contribution might be
negligible in the far-infrared wavelengths like 71.42 (hereafter 70) $\mu$m by
Spitzer (Lellouch et al., 2011) or longer wavelengths by Herschel (Lellouch et
al., 2016). This effect is important if the haze particles behave optically
similar to Titan-like tholins (Zhang et al., 2017). Recent studies suggested
that the haze particles might be primarily composed of hydrocarbon and nitrile
ices, which might produce a much smaller cooling effect (Lavvas et al., 2021).
A later study found that the ice particles could also explain the cold
temperature of Pluto’s atmosphere (Wan et al., 2021). In this case, the
infrared flux from the icy haze seems much smaller than the Titan-like organic
haze, but its contribution to the thermal lightcurve might not be negligible.
The large thermal emission from haze could obscure the surface features and
greatly reduce the relative amplitude of rotational lightcurve variations in
the mid-infrared. As a result, including the haze emission in the rotational
lightcurve analysis would not only lead to a better understanding of how the
haze affects the derivation of surface properties from observed lightcurves,
but also provide an opportunity to constrain the haze emission, which is
largely uncertain in atmospheric models due to the lack of constraints from
composition and optical constants for icy particles (Lavvas et al. 2021,
Section 2.2 and Discussion in Wan et al. 2021) and non-icy tholin-like
materials (Khare et al., 1984; Jovanović et al., 2020, 2021; Moran et al.,
2022).

To better constrain both the surface properties and haze emission, in this
study we also introduce a Bayesian retrieval framework that can efficiently
evaluate the posterior distribution of parameters. The Bayesian inference
technique has been widely used in analyzing remote sensing data on solar
system planets (e.g., Zhang et al., 2013; Fan et al., 2019) and exoplanets
(e.g., Line et al., 2011, 2013). It has also been used to retrieve the
chemical abundances in Pluto’s atmosphere from occultation data (Young et al.,
2018) as well as the haze distribution like particle sizes (Fan et al., 2022).
However, it has not been applied to Pluto’s rotational lightcurve analysis
before. With detailed surface albedo maps of Pluto and Charon now available,
using Bayesian inference might provide useful constraints on the posterior
distribution of the derived surface properties and haze emission.

We analyze three rotational lightcurves of the Pluto-Charon system in the
infrared, selecting data from Spitzer at 24 and 70 $\mu$m in 2004 (Lellouch et
al., 2011), as well as data from Herschel at 70 $\mu$m in 2012 (Lellouch et
al., 2016). We select the 24 $\mu$m data as the haze emissions could be
significant at this wavelength (Zhang et al., 2017), while the 70 $\mu$m data
was selected as haze emissions are not important at this wavelength (as
explained in Section 2.3). By using both channels, we wish to differentiate
the contribution of haze in the lightcurves. We include the Herschel data
because of its better quality compared to the Spitzer data. However, we do not
include the Herschel data at sub-mm wavelengths (such as 500 $\mu$m, Lellouch
et al., 2016) due to the significant subsurface contribution and seasonal
effects that could complicate the haze flux retrieval analysis.

Our paper is organized as follows: In Section 2, we introduce the surface maps
of Pluto and Charon from New Horizons, our thermal physical model, as well as
the retrieval models; In Section 3, we discuss our retrieval results in two
scenarios, with and without the haze contribution. We also use our best-fit
models to predict the rotational lightcurves of several channels on JWST’s
Mid-Infrared Instrument (MIRI) in Section 3; Finally, we conclude this paper
in Section 4.

##  2 Methods

###  2.1 Maps of Albedo and Emissivity

Figure 1: Surface maps of Pluto and Charon derived from New Horizons
observations. Part of the southern hemispheres were in winter darkness (polar
night) during the flyby, and thus the albedo distributions are not well
determined south of $\sim$ -33∘. (a) Map of Pluto’s bolometric hemispherical
albedo from Hofgartner et al. (2023). (b) Histogram of Pluto’s albedo
distribution, where local minima of 0.36 and 0.82 are adopted as unit
boundaries. (c) The distribution of different albedo units on Pluto’s surface.
(d) Map of Charon’s albedo from Buratti et al. (2017). Note that, Charon
rotates synchronously with Pluto and their lit hemispheres exhibit a 180∘
longitude difference.

The surface temperatures of Pluto and Charon are greatly affected by albedo,
emissivity, and thermal inertia. The distributions of all three parameters on
Pluto are primarily determined by the distribution of volatile ices and non-
volatile materials (Grundy et al., 2016; Moore et al., 2016; Schmitt et al.,
2017). Surfaces of the Pluto-Charon system were observed during the New
Horizons flyby in 2015, including high-resolution mosaics of the “near side”
and $\sim$20$\times$ poorer imaging of the “far side” (Buratti et al., 2017;
Stern et al., 2021).

As a disk-integrated quantity, Bond albedo was widely used to measure the
global energy balance of a planetary body. To calculate the local surface
temperature, here we use the bolometric, hemispherical (disk-resolved) albedo.
Preliminary maps of disk-resolved albedo for Pluto and Charon using
panchromatic LORRI observations were presented in Buratti et al. (2017).
Hofgartner et al. (2023) updated the Pluto albedo map by fitting the lunar-
Lambert photometric function (e.g., Buratti & Veverka, 1983) to both LORRI and
MVIC measurements from a wide range of imaging geometries and by considering
different photometric functions for Pluto’s extreme dark and extreme bright
terrains. The bolometric albedo was approximated using New Horizons LORRI
panchromatic and MVIC panchromatic and color filter observations. The updated
albedo map of Pluto is shown in Figure 1(a) and the map of Charon is shown in
Figure 1(d).

In our study, we directly apply these latest albedo maps of Pluto and Charon
from the New Horizons 2015 flyby to explain the 2004 Spitzer data and 2012
Herschel data. In other words, we assume the albedo distribution does not
change between 2004 and 2015. For Pluto, this is only an approximation because
the ice distribution on its surface might change with season (Bertrand &
Forget, 2016; Bertrand et al., 2020). However, this is the only high-
resolution map of Pluto’s surface to date. Before the New Horizons flyby,
Pluto’s surface was barely resolved. Thus, Lellouch et al. (2011) chose to use
Pluto’s HST reflection lightcurve in 2002–2003 to constrain the albedo of
different geological units in their assumed surface maps. One might achieve a
more accurate albedo map of Pluto in 2003 based on both the New Horizons and
the HST data. That is beyond the scope of this study as it requires a detailed
investigation using an atmosphere-surface model of Pluto (e.g., Young, 2017;
Bertrand et al., 2020). As shown in Section 3, using the New Horizons albedo
maps can produce a decent fit of both the 2004 Spitzer and 2012 Herschel
observations in our retrieval.

Due to the observing geometry, the New Horizons maps do not cover locations
south of $\sim$ -33∘ for both Pluto and Charon. Because the sub-solar and sub-
observer latitudes were 34.5∘ and 32.9∘ in 2004 (Spitzer) and about 47.0∘ in
2012 (Herschel), the missing albedos in the southern hemisphere would not
strongly impact the lightcurves due to the projection effect. We tested a
range of albedo values from 0 to 1 in those regions, and confirmed that our
results were not changed. In our study, we utilize 120$\times$60 congrid maps
for both Pluto and Charon, corresponding to a spatial resolution of 3∘. We
also tested lightcurve simulations with higher spatial resolution, and the
results were unchanged.

The emissivity distributions of Pluto and Charon are not directly constrained
by New Horizons observations. For simplicity, we follow the approach in
Lellouch et al. (2011) that separated Pluto’s surface into three geological
units based on albedo and assumed a constant emissivity for each unit. We
associate such composition units: dark H2O ice/tholin mix (unit 3), CH4 ice of
intermediate albedo (unit 2), and bright N2 ice (unit 1) with three albedo
units in Figure 1(b) based on the histogram of Pluto’s updated albedo
distribution from New Horizons: unit 3 of darkest ices (0$-$0.36), unit 2 of
intermediate ices (0.36$-$0.82), and unit 1 of brightest ices (0.82$-$1.0).
The divided units (Figure 1c) appear to represent an approximate distribution
of icy components identified by LEISA (Grundy et al., 2016; Schmitt et al.,
2017; Gabasova et al., 2021; Emran et al., 2023), although multiple ice
mixtures show a much more complex pattern in the infrared LEISA maps. In
contrast, Charon’s surface (Figure 1d) is relatively uniform, and we use only
one emissivity value. The thermal inertia distributions are also less
constrained than that of albedo and we use a single value for Pluto and a
separate single value for Charon.

###  2.2 Thermophysical Model of Surface

Similar to previous work (Spencer et al., 1989; Lellouch et al., 2000, 2011;
Forget et al., 2017), we model the surface and subsurface temperature
distribution considering incoming solar heating, surface thermal emission, and
thermal conduction for Pluto and Charon. Sublimation and deposition of
volatiles are not included. At every point on the surface, the temperature
($T$) of subsurface layers at a time ($t$) and depth ($z$) obeys a diffusion
equation:

$\rho c_{p}\frac{\partial T}{\partial t}=\frac{\partial}{\partial
z}({k\frac{\partial T}{\partial z}}),$ (1)

where $\rho$ (kg m-3) is the density, $c_{p}$ (J K-1 kg-1) is the specific
heat capacity and $k$ (W m-1 K-1) is the thermal conductivity. $\rho c_{p}$ =
106 J m-3 K-1 is adopted (Forget et al., 2017) everywhere at any depth on
Pluto. The upper ($z=0$) and lower ($z=\infty$) boundary conditions are

$\displaystyle k\frac{\partial T}{\partial z}|_{z=0}+Q_{s}=E_{s},$ (2)
$\displaystyle\frac{\partial T}{\partial z}|_{z=\infty}=0,$ (3)

respectively. At the surface, thermal conduction and the absorbed insolation
($Q_{s}$) balance with the infrared emission to space ($E_{s}$) (Hayne et al.,
2017). We assumed no heat flux at the lower boundary.

To intuitively understand how deep the surface heat flux can propagate
downwards to affect the subsurface temperature, we can introduce the skin
depth ($l_{s}$)

$l_{s}=\frac{\Gamma}{\rho c_{p}}\sqrt{\frac{P}{\pi}},$ (4)

where $P$ is the period of stellar flux forcing and $\Gamma=\sqrt{k\rho
c_{p}}$ is the thermal inertia. The depth of 8$l_{s}$ approximates the lower
boundary, below which the periodic insolation changes hardly impact the
substrate temperature (Spencer et al., 1989).

Pluto has two dominant periods. For the diurnal cycle, $P_{day}$ = 6.4 Earth-
days. If we assume the thermal inertia in the shallow subsurface is
$\Gamma_{day}$ = 20 MKS (J s${}^{-\frac{1}{2}}$ m-2 K-1) from Forget et al.
(2017), the corresponding thermal diffusivity is $k$ = 4$\times$10-4 W m-1 K-1
and the lower boundary depth is $8l_{s}$ = 0.067 m. For the seasonal cycle,
$P_{year}$ = 248 Earth-years. Thermal inertia of the deep subsurface
$\Gamma_{year}$ = 800 MKS (Forget et al., 2017) corresponds to $k$ = 0.64 W
m-1 K-1 and $8l_{s}$ = 320 m, much deeper than the diurnal skin depth. Pluto
is assumed to have homogeneous thermal inertia and conductivity; these values
of $\Gamma$ and $k$ are used at all locations. Because our study focuses on
rotational lightcurves, we only consider the shallow subsurface for the
diurnal cycle in a self-rotating model, in line with Spencer et al. (1989). In
Appendix A, we also test both diurnal and seasonal forcing using a two-layer,
orbit-rotating model, which has a shallow subsurface with low thermal inertia
overlaying a deep subsurface with high thermal inertia. The full-orbit
simulations using this two-layer model show that the diurnal surface
temperature variation in our self-rotating model is not greatly impacted by
the seasonal cycle.

Evolution of the temperature is solved with a central difference explicit
method and can be discretized by

$T_{j}^{n+1}=T_{j}^{n}+dt\frac{k_{j}}{(\rho
c_{p})_{j}}\frac{1}{(z_{j+1}-z_{j-1})/2}(\frac{T_{j+1}^{n}-T_{j}^{n}}{z_{j+1}-z_{j}}+\frac{T_{j}^{n}-T_{j-1}^{n}}{z_{j}-z_{j-1}}),$
(5)

where the subscripts of $j$ and $n$ represent different levels and numerical
time steps, respectively. Similar to Hayne et al. (2017), the boundary
conditions with depth can be discretized by

$\displaystyle
k_{0}\frac{-3T_{0}^{n+1}+4T_{1}^{n+1}-T_{2}^{n+1}}{2(z_{1}-z_{0})}=\epsilon\sigma{T_{0}^{n+1}}^{4}-(1-A)F_{0}^{n+1},$
(6) $\displaystyle T_{N}^{n+1}=T_{N-1}^{n+1},$ (7)

where $A$ is the albedo (bolometric hemispherical albedo for Pluto and disk-
resolved approximate Bond albedo for Charon), $\epsilon$ is the bolometric
emissivity, $\sigma$ is the Stefan-Boltzmann constant, and $F_{0}$ (W m-2) is
the insolation that varies with time and location. For a given time and
orbital phase, we first calculate the sub-solar point on Pluto using Equation
(15) in Ohno & Zhang (2019) and then the value of $F_{0}$ at any surface point
based on its incident angle.

The model is divided into 15 uneven vertical layers, with the first two
starting at a thickness of 10-3 m, below which the thickness increases by a
factor of 1.5 between adjacent layers. The total depth of the model is around
0.3 m, much deeper than $8l_{s}$ (less than 0.1 m) for $\Gamma$ ranging from
20 to 30, as suggested in Lellouch et al. (2011). The thermal conduction
timescale near the surface is $\sim$625 s, increasing with depth due to larger
layer thickness. To accurately resolve a Pluto day during rotation, we set
1800 time steps with a time step of $\sim$300 s. The initial temperature
profile is isothermal based on the surface’s equilibrium temperature under
diurnal-mean insolation. For each point on the surface, we run the model under
diurnal insolation forcing until the diurnal variation of the surface
temperature does not change between two subsequent rotations. Typically, five
rotations are sufficient to reach convergence.

Given the albedo maps, emissivity, and thermal inertia, we can calculate the
temporal and spatial distribution of the surface temperature of Pluto and
Charon with our thermophysical model. However, the model does not consider the
ice sublimation-condensation processes, which might impact the surface ice
temperature through the latent heat exchange. Following Lellouch et al.
(2011), we include additional treatment of N2 and CH4 ice temperature in our
model. For the N2 ice, given its high albedo (Figure 1) and bolometric
emissivity of 0.5 (Lellouch et al., 2011), the N2 ice radiative-equilibrium
temperature is low and is mainly controlled by the sublimation-condensation
process. We can just fix the N2 ice temperature to its vapor-ice equilibrium
in our model. Based on the surface pressure of Pluto in 2004 and 2012 from
stellar occultation and modeling (Meza et al., 2019) and the saturation vapor
pressure of the N2 ice in the beta form (Fray & Schmitt, 2009), we estimate
the N2 ice temperature to be 36.35 K in 2004 and 36.86 K in 2012. For the CH4
ice, it was found that the sublimation cooling effect becomes important when
the CH4 ice temperature exceeds $\sim$54 K (Stansberry et al., 1996b). Below
54 K, the radiative process controls the surface temperature. Thus, we set a
maximum temperature of 54 K for the CH4 ice, which was also adopted in
Lellouch et al. (2000, 2011).

###  2.3 Outgoing Thermal Emission

To model the rotational lightcurves, we first define the sub-observer point
($\alpha_{0}$, $\delta_{0}$), where $\delta_{0}$ is the latitude and
$\alpha_{0}=\omega t$ is the longitude at time $t$ with a rotation rate
$\omega=2\pi/P_{day}$. The outgoing surface emission ($F_{s}$) from Pluto and
Charon at a given wavelength ($\lambda$), when viewed from the sub-observer
point, can be calculated as follows:

$F_{s,\lambda}(\alpha_{0}(t),\delta_{0})=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}\epsilon_{\lambda}(\alpha,\delta)B_{\lambda,T}(\alpha,\delta)\cos(\delta)\cos(\Theta)d\delta
d\alpha,$ (8)

where $\epsilon_{\lambda}$ is the spectral thermal emissivity at wavelength
$\lambda$, and $B$ (W m-2 Sr-1 $\mu$m-1) is the blackbody emission determined
by the surface temperature at a surface location ($\alpha$, $\delta$).
$\cos(\Theta)=\max[\sin(\delta_{0})\sin(\delta)+\cos(\delta_{0})\cos(\delta)\cos(\alpha-\alpha_{0}),0]$,
where $\Theta$ means the angular distance between the sub-observer point
($\alpha_{0}$, $\delta_{0}$) and surface locations ($\alpha$, $\delta$). We
generate the rotational emission lightcurves by using a fixed sub-observer
latitude ($\delta_{0}$) and varying the time $t$ from 0 to $2{\pi}/{\omega}$.

Previous thermal emission calculations only considered the contribution from
Pluto’s and Charon’s surfaces (Lellouch et al., 2011) but neglected Pluto’s
haze flux, which might significantly affect the total outgoing emission.
Pluto’s haze extends several hundred kilometers above the surface, and the
particles can scatter and absorb solar radiation in the optical wavelengths
and emit in the mid-infrared. Due to the uncertainty in haze composition, the
estimated haze infrared flux at 24 $\mu$m from atmospheric models differs by
more than an order of magnitude, ranging from 0.1$\times$ to a few mJy (Zhang
et al., 2017; Lavvas et al., 2021; Wan et al., 2021). The latter is comparable
to the surface emission of Pluto at the same wavelength (Lellouch et al.,
2011). Thus, it is necessary to include Pluto’s haze when analyzing the
infrared spectra and thermal lightcurves of the Pluto-Charon system.

We need to consider several facts to account for the effects of Pluto’s haze
in the emission calculation. First, the haze might block part of the thermal
emission from the surface, obscuring the surface features. However, the
infrared optical depth of the haze is too small ($\sim 10^{-3}$ at 24 $\mu$m,
Zhang et al. 2017) to be important. Therefore, we can directly add the haze
emission to the surface emission in our model.

Second, the haze emission contributes differently to the total outgoing flux
at different wavelengths. The haze emission might be as important as the
surface emission in the mid-infrared like 24 $\mu$m, but it is negligible in
the far-infrared like 70 $\mu$m. The haze flux, calculated assuming Titan-like
tholins in Zhang et al. (2017), is about 2 mJy at 70 $\mu$m, much smaller than
the blackbody emission of Pluto’s surface at this wavelength, which is over
100 mJy. As a result, in this study, we only consider the haze emission at
shorter wavelengths such as 18 or 24 $\mu$m but not at longer wavelengths such
as 70 $\mu$m.

Lastly, Pluto’s global haze layers are not entirely spatially uniform (Cheng
et al., 2017) as shown in the north-south brightness asymmetry in the LORRI
images. 3D haze transport simulation on Pluto (Bertrand & Forget, 2017) also
shows that more haze was seen at the North pole due to the larger incoming
solar flux. However, Pluto’s gases and haze are predicted to be quickly mixed
over all longitudes by the zonal circulation, therefore justifying the
approximation of a longitudinally uniform haze layer. Because this study
focuses on the disk-averaged flux, we assume a constant haze emission when
Pluto is rotating. The observed inverted temperature profile (Hinson et al.,
2017) above Pluto’s cold icy surface was classically understood as being
primarily caused by thermal conduction and CH4 heating/cooling radiative
equilibrium, with some additional role due to CO, HCN, and C2H2 cooling (Yelle
& Lunine, 1989; Strobel et al., 1996; Strobel & Zhu, 2017). However, it has
recently been shown that haze radiative heating/cooling (Zhang et al., 2017)
and eddy heat transport (Wan et al., 2021) could dominate the region. If haze
is the main coolant, the local surface temperature does not strongly influence
the haze emission whose contribution mainly comes from 20 km above the surface
where the atmospheric temperature maximizes (Hinson et al., 2017).

###  2.4 Bayesian Inference Retrieval

Directly solving for the surface properties of Pluto and Charon, as well as
Pluto’s haze flux, is a challenging task due to the complexity of the
thermophysical model and lightcurve calculations involved. We choose a
Bayesian retrieval method, which has been widely used in planetary and
exoplanet studies, to simultaneously constrain all parameters and assess their
uncertainties. This method iteratively solves the inverse problem to find the
best-fit model solution to the observed data. Various algorithms have been
developed to fit observational data and determine the model parameters as well
as their posterior distributions, such as Markov Chain Monte Carlo (Mackay,
2003) and Nested Sampling technique (Skilling, 2004). For an efficient and
unbiased sampling of the parameter space, we employ the well-developed
software package PyMultiNest (Buchner et al., 2014), which utilizes the Nested
Sampling technique. By comparing the simulated flux with observations from
Spitzer and Herschel in Lellouch et al. (2011, 2016), we are able to find the
best-fit model parameters that minimize the difference between the simulated
and observed lightcurves. In addition, the Nested Sampling approach can
estimate the posterior distribution (i.e., uncertainty) of the retrieved
parameters.

We list important input parameters of our model in Table 1, including the
distance and locations of the Sun and observer to Pluto in 2004 and 2012
(Lellouch et al., 2011, 2016), the radii of Pluto and Charon (Nimmo et al.,
2017), and the thermal properties of ices. The retrieved parameters in our
model include Pluto’s thermal inertia ($\Gamma_{Pluto}$), emissivity of unit 2
(approximately CH4 ice, $\epsilon_{CH_{4}}$), Charon’s thermal inertia
($\Gamma_{Charon}$), as well as the haze flux at 24 $\mu$m ($F_{haze}$). We
neglect the haze emission contribution at 70 $\mu$m as it is much smaller
compared to the surface flux at that wavelength. Because the N2 ice
temperature is low, its contribution to the rotational lightcurves is
negligible. Thus, we fix its bolometric and spectral emissivity using the
values in Lellouch et al. (2011). On the other hand, Pluto’s H2O ice/tholin
mix unit (dark terrain) and Charon’s surface generally have a high emissivity
Lellouch et al. (2011), so we fix their emissivity as unity in our model. We
assume the spectral emissivity of the CH4 ice at 24 and 70 $\mu$m is the same
as the bolometric emissivity. In other words, there is only one emissivity
value for the CH4 ice. It appears that one emissivity parameter of the CH4 ice
unit in our retrieval model is sufficient to explain the observed lightcurves.
Even though the CH4 emissivity might be wavelength-dependent (Stansberry et
al., 1996a) and previous studies (Lellouch et al., 2011) tested wavelength-
dependent emissivity, the current data could not constrain all three
independent values.

| Spitzer (2004) | Herschel (2012) | JWST (2023)
---|---|---|---
Sub-solar distance (AU) | 30.847 | 32.19 | 34.6
Sub-solar latitude (∘) | 34.5 | 47.0 | 58.3
Sub-observer distance (AU) | 30.95 | 32.4 | 34.4
Sub-observer latitude (∘) | 32.9 | 48.5 | 57.3
Planetary radius (km) | 1188 for Pluto and 606 for Charon
Emissivity for N2 ice | 0.5 for bolometric and 70 $\mu$m, 0.05 for others∗
Emissivity for H2O ice/tholin mix | 1.0 for bolometric and all wavelengths
Emissivity for Charon | 1.0 for bolometric and all wavelengths
Fixed temperature for N2 ice (K)† | 36.35 | 36.86 | 36.35
Max temperature for CH4 ice (K) | 54
Table 1: Input parameters of the thermophysical model. ∗Based on the N2
emissivities in Lellouch et al. (2011), we adopted 0.5 for bolometric
emissivity and the emissivity at 70 $\mu$m, and 0.05 for other shorter
wavelengths between 18 and 30 $\mu$m. †Due to the very low temperature and
small spectral emissivity, this value does not affect Pluto’s total thermal
emission.

In the Nested Sampling approach, live points are the set of samples being
actively updated by the algorithm during the course of the computation. At
each iteration of sampling, the point with the lowest likelihood value (i.e.,
the least probable point) is removed from the set of live points and replaced
with a new point that has a higher likelihood value. This process continues
until the desired level of precision is achieved or until the algorithm
converges. Increasing the number of live points can lead to more accurate
estimates of parameters and their posterior distributions, but also increases
the computational cost. The number of live points should be at least 10 times
and preferably 25 times the number of free parameters. We choose 1600 live
points for our retrieval, although 800 also works and gives similar results.

Figure 2: The posterior distribution of retrieved parameters for the “free
retrieval” scenario with haze (red) and “surface only” scenario without haze
(cyan), including 24 $\mu$m haze flux $F_{haze}$ (mJy), CH4 ice emissivity
$\epsilon$, thermal inertia $\Gamma$ (MKS) for both Pluto and Charon. Vertical
dashed lines in the diagonal panels denote 2.5%, 50%, and 97.5% quantiles of
the distributions. Off-diagonal panels show the two-dimensional parameter
covariance plots, indicating a strong negative correlation between the haze
flux and Charon’s thermal inertia (lower-left panel). Figure 3: Rotational
emission lightcurves of Pluto and Charon from our best-fit models (lines) and
observations (dots) at Spitzer 24 $\mu$m and 70 $\mu$m (first two columns),
Herschel 70 $\mu$m (third column). Contributions from Pluto’s different
components are separately plotted for the “free retrieval” scenario with haze
(upper panel) and “surface only” scenario (lower panel). Figure 4: Surface
temperature distributions of Pluto (left) and Charon (right) in 2004 from our
best-fit models in the “free retrieval” scenario with haze (upper panel) and
the “surface only” scenario (lower panel). We plotted the maps at noon time
when the sub-solar point at Pluto is 180∘ East Longitude. The animations of
diurnal temperature evolution in 2004 and 2012 are available on the journal’s
website.

##  3 Results

###  3.1 Two Possible Scenarios to Explain Lightcurves

We first present the results from the “free” retrieval, where we include the
haze flux as a free parameter to constrain its value based on the observed
lightcurves. The resulting posterior distribution of retrieved parameters and
the fitted rotational lightcurves are presented in Figure 2 (red) and Figure 3
(upper panel), respectively.

Using the best-fit parameters, our model is able to explain the rotational
lightcurves relatively well at both 24 and 70 $\mu$m. Both Pluto and Charon
emit significantly more at 70 $\mu$m than at 24 $\mu$m, as shown in Figure 3.
The corresponding temperature exhibits significant diurnal and spatial
variation. We show the snapshots of their temperature maps at the noon time in
Figure 4. We also provide animations of the diurnally-varying temperature
distributions of Pluto and Charon from our best-fit models, which are
available on the journal’s website.

The lightcurves (black in Figure 3) at 24 and 70 $\mu$m exhibit two peaks,
with the main peak occurring near 100∘ and a smaller peak near 300∘ in
longitude. The main peak is primarily determined by the distribution of the
dark and warm H2O ice/tholin mix (purple) on Pluto, which is unevenly
distributed along the longitude at low latitude regions (unit 3 in Figure 1).
One can also see the corresponding temperature map at noon time when during
Spitzer’s observations in 2004 (Figure 4a). While the CH4 ice, which covers a
larger surface area in the northern latitudes (unit 2 in Figure 1), emits
comparable flux to that of the H2O ice/tholin mix but shows less variation
across longitude. This CH4 ice unit contributes to the second but smaller
emission peak near 300∘ in longitude. In contrast, the N2 ice emission (green)
is negligible due to its much lower temperature. As the lightcurve variation
primarily arises from Pluto’s surface, it provides strong constraints on
parameters such as Pluto’s thermal inertia and emissivity. Unlike Pluto,
Charon’s uniform but darker surface (blue) produces almost flat lightcurves of
about 2.2 mJy at 24 $\mu$m and 128.5 mJy at 70 $\mu$m, which is smaller than
Pluto’s total surface emission (gray).

In the best-fit model, Pluto’s haze that was ignored in previous work
(Lellouch et al., 2000, 2011, 2016), emits comparably to Charon’s emission at
24 $\mu$m. The haze emission of 1.63 mJy at 24 $\mu$m is similar to that in
the atmospheric model calculations using Titan-like tholins in Zhang et al.
(2017). However, the posterior distribution of this parameter indicates that
the retrieved haze flux is 1.63${}^{+0.43}_{-0.49}$ mJy, given by live points
distributed within the 2.5% and 97.5% quantiles of the posterior by the
Nesting Sampling technique (Figure 2). The lower value of 1.14 mJy at 24
$\mu$m might correspond to the situation of less absorbing particles in Zhang
et al. (2017). Note that, in our model the haze contribution is negligible at
70 $\mu$m (see Section 2.3).

The wide range of retrieved haze flux on Pluto is due to the inherent
degeneracy between Pluto’s haze emission and Charon’s surface emission. This
is because Spitzer and Herschel telescopes cannot resolve Charon from Pluto
and both Pluto’s haze and Charon’s surface are spatially uniform. Even though
we have tried to include the haze emission at 24 $\mu$m but not at 70 $\mu$m,
and Charon contributes to both wavelengths, the current data is insufficient
to break the degeneracy. As illustrated in the corner plot that shows the
joint distribution of the parameters (lower-left panel of Figure 2), Pluto’s
haze emission exhibits a strong negative correlation with Charon’s thermal
inertia. For instance, when the haze flux is approximately 1.14 mJy at 24
$\mu$m, the retrieved Charon’s thermal inertia is about 13 MKS, yielding about
2.87 mJy. If the haze flux is as high as 2.06 mJy, the corresponding derived
Charon’s thermal inertia is 45 MKS, corresponding to a Charon surface flux of
about 1.75 mJy. This correlation implies that the total emission from Pluto’s
haze and Charon’s surface is roughly constant in all solutions, which is about
4 mJy at 24 $\mu$m (Figure 3).

Moreover, in any retrieval process, the goodness of fit is influenced not only
by data quality but also by the accuracy of the forward model. Our
thermophysical model has made several assumptions, including dividing Pluto’s
surface into three geological units and a single thermal inertia value for the
entire surface. Due to these model assumptions and incomplete information
regarding the input albedo distribution (and neglecting its temporal
evolution), our best-fit model does not perfectly match all data points in the
observations, particularly those with minimal error bars (e.g., Figures 3b and
3e). Recognizing the potential biases in the forward model’s construction and
the substantial uncertainty range in the retrieved haze emission flux, we
should be careful when interpreting the posterior distribution of the
retrieved parameters. As an extreme scenario, it would be interesting to see
the best fit using a model without haze. Therefore, we conduct an additional
retrieval in a “surface-only” scenario, where the haze flux is excluded, and
only surface properties are retrieved. This approach resembles the one in
Lellouch et al. (2011), albeit with updated albedo maps from New Horizons.

The parameters and their posterior distributions from the “surface only”
retrieval are shown in Figure 2 (cyan). The fitted rotational lightcurves are
shown in Figure 3 (lower panel). Interestingly, our results are consistent
with that in Lellouch et al. (2011). The “surface only” scenario using the New
Horizons albedo maps can also fit the rotational lightcurves relatively well
without any haze contribution. As expected, the “surface only” retrieval
yields low thermal inertia of Charon (5.5$\pm$0.5 MKS), much lower than the
“free retrieval” results and previous estimates when using a different albedo
map of Pluto (Lellouch et al., 2011, 2016). It is still acceptable because
trans-Neptunian objects (TNOs) have even smaller values of 2.5$\pm$0.5 MKS
(Lellouch et al., 2013). The contribution of Charon’s emission flux at 24
$\mu$m is about 4.2 mJy in the “surface only” case, roughly equal to the total
emission of both Pluto’s haze and Charon’s surface in the “free retrieval”
case. Given that the “free retrieval” scenario prefers a non-zero haze flux
and that the “surface only” scenario can also provide a decent fit, we
conclude that the current data cannot constrain Pluto’s haze flux. It also
means that the current data cannot distinguish whether Pluto’s haze is Titan-
like tholins (Zhang et al., 2017) or ice particles (Lavvas et al., 2021; Wan
et al., 2021).

Radio observations of the Pluto-Charon system have a high spatial resolution
that allows for the separation of Charon from Pluto. For example, SMA observed
the dayside brightness temperature of Charon at 1.36 mm, yielding a
temperature of 43$\pm$13 K in 2005. Our best-fit models give Charon’s surface
brightness temperature of 51.1 K for the “free retrieval” scenario with haze
and 54.0 K for the “surface only” scenario without haze, respectively. Both
temperatures are within the uncertain range of the SMA results. ALMA
observations at 0.86 mm showed a temperature of 43.9$\pm$0.3 K for Charon in
2015, whereas our models predict a temperature of 51.8 K for the “free
retrieval” scenario with haze and 54.3 K for the “surface only” scenario
without haze. These derived temperatures are greater in the scenario absent of
haze due to smaller thermal inertia from retrievals, however, temperatures
from both scenarios are larger than the corresponding radio observations.
Similar results are also found for other observations, such as SMA at 1.10 mm
in 2010 and VLA at 9.00 mm in 2011 (see Table 10 in Bird et al. 2019), New
Horizons dayside-staring at 4.17 cm in 2015 (Bird et al., 2019), and Herschel
at 500 $\mu$m in 2012 (Lellouch et al., 2016). This discrepancy is expected
because radio observations do not directly probe the surface temperature but
rather a mix of contributions from the surface and subsurface, which could be
colder than the surface on the dayside. Furthermore, the radio emissivity of
surfaces is not unity (Lellouch et al., 2016). Therefore, current radio
observations of Charon’s surface cannot distinguish the “free retrieval”
scenario with haze versus the “surface only” scenario.

Observations of near-infrared spectra can also be used to probe the surface
ice temperature, such as using the absorption bands of H2O ice as thermometers
(Grundy et al., 1999; Protopapa et al., 2021). Buie & Grundy (2000) inferred a
disk-average temperature of 60$\pm$20 K for Charon’s surface via this method.
Recently, Holler et al. (2017) revisited this effort on the 1.65 $\mu$m
absorption band of H2O ice and reported Charon’s mean surface temperature as
45$\pm$14 K in 2015. Our modeled surface temperature for both scenarios (49.5
K with haze and 48.9 K without haze in 2012) falls within this large
uncertainty range, implying that the current near-infrared “thermometer” is
still too large to precisely constrain Charon’s surface temperature to
separate the two scenarios in this study.

While we cannot well constrain the properties of Pluto’s haze and Charon’s
surface, we can determine their total emission flux at 24 $\mu$m in 2004,
which is about 4 mJy in both haze and haze-free scenarios. In other words, the
emission from Pluto’s surface can be constrained. Indeed, Pluto’s surface
properties appear well-determined from the shape of lightcurves. The surface
temperature distributions of Pluto in the best-fit models appear very similar
between the “free retrieval” scenario with haze (Figure 4a) and the “surface
only” scenario without haze (Figure 4c). In both scenarios, the retrieved
bolometric emissivity of Pluto’s CH4 ice is about 0.4, suggesting that the CH4
ice surface might have small grains according to Figure 5 in Stansberry et al.
(1996a). Our retrieved CH4 emissivity is smaller than that from Lellouch et
al. (2011, 2016), which yields the bolometric emissivity of 0.7–0.9, spectral
emissivity of 1.0 at 24 $\mu$m and 0.9 at 70 $\mu$m. The derived thermal
inertia of Pluto is 10–20 MKS in Lellouch et al. (2011) and 16–26 MKS in
Lellouch et al. (2016), higher than our optimal value (8–10 MKS). The reason
has yet to be explored but is likely related to the difference in albedo maps
in the two studies.

In summary, incorporating haze flux into retrievals is consistent with all
existing observations. However, the current data quality is insufficient to
distinguish between the contributions of Charon and haze, and cannot directly
constrain the haze emission of Pluto. In the next section, we suggest that
future observations using the JWST at additional wavelengths may assist in
resolving this issue.

Figure 5: Rotational emission lightcurves of Pluto and Charon simulated for
the JWST observations at 18 $\mu$m, 21 $\mu$m, and 25.5 $\mu$m. Results of the
“free retrieval” scenario with haze (upper panel) and the “surface only”
scenario without haze (lower panel) are plotted with contributions from
different geologic units and haze on Pluto and the surface of Charon. The haze
emission fluxes at multiple wavelengths are calculated using the atmospheric
model from Zhang et al. (2017) assuming Titan-like tholins. Figure 6:
Variation of rotational emission lightcurves (in units of $\mu$Jy) simulated
for JWST’s observations at 18 $\mu$m, 21 $\mu$m, and 25.5 $\mu$m. Results of
the “free retrieval” scenario with haze (upper panel) and the “surface only”
scenario without haze (lower panel) are plotted. Dashed lines indicate the
10$\sigma$ faint-source detection limits in a 104-second integration for
corresponding channels of the JWST/MIRI.

###  3.2 Predictions for JWST

Future observations with JWST can help distinguish between the two scenarios
and provide valuable information on the haze and surface properties of Pluto.
To achieve this, we have selected three typical channels of the JWST/MIRI: 18,
21, and 25.5 $\mu$m, with sensitivities (Glasse et al., 2015) of 4.87, 9.15,
and 30.8 $\mu$Jy, respectively. These channels are chosen as we primarily
focus on the thermal emission from the Pluto-Charon system, and shorter
wavelengths may be contaminated by reflected light. To support the planned
JWST observations of Pluto in 2022–2023 111https://www.stsci.edu/jwst/science-
execution/program-information.html?id=1658, we calculate the rotational
lightcurves based on a sub-solar latitude of 58.3∘, a sub-observer latitude of
57.3∘, a distance of 34.6 AU to the Sun, and a distance of 34.4 AU to the
observer of JWST. In this configuration, it is possible to resolve Charon from
Pluto and break the degeneracy between Pluto’s haze emission and Charon’s
surface emission. Two limiting cases without haze and with thick haze
(assuming Titan-like tholins) are presented in Figure 5 using the best-fit
parameters in the “surface only” and “free retrieval” scenarios in Section
3.1, respectively. We assume the spectral emissivity of the surface units at
these three wavelengths is the same as that at 24 $\mu$m from the retrieval.

In the “surface only” scenario (lower panel in Figure 5), the surface emission
from both Pluto and Charon increases with wavelength as it approaches the
blackbody peak. At 18 $\mu$m, Pluto’s surface emission is 0.1 mJy, which
increases to 0.5 mJy at 21 $\mu$m and to 3.0 mJy at 25.5 $\mu$m. Similarly,
Charon’s surface emission increases from 0.25 mJy at 18 $\mu$m to 1.0 mJy at
21 $\mu$m and to around 4.5 mJy at 25.5 $\mu$m. Notably, Charon’s surface
emission is approximately twice larger than that of Pluto’s surface for all
three wavelengths. The contributions of different surface ices to Pluto’s
surface emission look dramatically different from those in 2004. This is
because the sub-solar and sub-observer latitudes have moved much further
north, from $\sim$34∘ in 2004 to $\sim$58∘ in 2023. Since most of the H2O
ice/tholin mix is concentrated near the equator and Pluto’s northern
hemisphere is predominantly covered with CH4 ice, the emission of Pluto’s
surface (gray) mainly comes from the CH4 ice (orange) rather than the H2O
ice/tholin mix (purple) in 2023 (see Figure 5). Therefore, the total
lightcurves (black) are less variable compared to that in 2004 and 2012
(Figure 3) and peak at 300∘ in 2023 rather than 100∘ in longitude.

Taking Pluto’s haze into consideration would make it more complicated to
interpret the data, but it also offers an opportunity to constrain Pluto’s
haze emission. If the haze is made of absorbing Titan-like tholins,
calculations using atmospheric model (Zhang et al., 2017) predict that it
could emit much more than Pluto’s surface at 18 and 21 $\mu$m (upper panel in
Figure 5). JWST observations at these two wavelengths can be used to directly
constrain the haze emission. On the other hand, if the haze emission is
smaller, the surface emission and haze emission cannot be easily disentangled.
As an extreme example, if the haze is made of brighter ice particles as
suggested by the chemical models (Lavvas et al., 2021), the haze radiative
cooling could be ten times smaller than that from the Titan-like tholins (Wan
et al., 2021). For instance, if the tholin haze opacity is reduced by 20 times
in the lower atmosphere as suggested by Wan et al. (2021), such haze emission
in the mid-infrared is comparable to the surface emission at 18 $\mu$m (0.1
mJy), a few times smaller than the surface at 21 $\mu$m, but negligible
compared with the surface at 25.5 $\mu$m. Thus it is best to use 18 $\mu$m to
constrain the haze properties. Yet, the emission from icy haze could also be
large within its solid-state absorption band, for instance, within the 10-15
micron range for nitrile ices (see Figure 4 in Zhang et al. 2017). Future
spectra from the JWST may be capable of identifying these spectral structures.
In any case, the rotational lightcurves from JWST in 2023 appear much flatter
than the Spitzer observations as sub-observer latitude moves further north in
absence of contribution from the H2O ice/tholin mix. It might be difficult to
distinguish the CH4 ice emission at high latitudes and the haze emission.
Combining with previous data from Spitzer or Herschel might be important to
break the degeneracy, although assumptions need to be made such as that the
ice distributions and haze temperature do not change significantly with time
from 2004 to 2023.

Interestingly, although the lightcurves of both Pluto and Charon exhibit
nearly flat profiles in 2023, their variations can be characterized using
JWST. In contrast to the lightcurves observed in 2004 and 2012, where Pluto’s
variation was much greater than that of Charon (Figure 3), our predictions
indicate that Pluto’s lightcurve variation during the JWST era will be
comparable to that of Charon (Figure 6). At wavelengths of 18, 21, and 25.5
$\mu$m, the peak-to-trough flux amplitudes are approximately 20, 100, and 400
$\mu$Jy, respectively (Figure 6), all of which exceed the JWST/MIRI
instrument’s detection limit. Pluto’s lightcurve exhibits a distinct “U”
shape, primarily influenced by the distribution of ice concentrated near the
180∘ mark in the northern hemisphere. On the other hand, Charon’s lightcurve
displays a maximum-flux phase shift from Pluto’s, likely influenced by
slightly darker materials in its northern high latitudes (Figure 1d) and the
thermal phase lag (note that both Pluto and Charon are retrograde rotators
with an obliquity larger than 90 degrees). Although distinguishing between
best-fit models with and without haze in the lightcurve variation may be
challenging (Figure 6), the shape of the lightcurve can serve as a useful
constraint for determining the thermal inertia and ice albedo distributions on
both bodies during the JWST era.

##  4 Summary and Discussions

In this paper, we investigate the properties of Pluto’s surface and haze using
thermophysical modeling and Bayesian retrieval methods. We re-analyze the
rotational lightcurves of the Pluto-Charon system observed by Spitzer and
Herschel, incorporating Pluto’s haze emission and the latest surface albedo
maps from New Horizons. Our main findings are as follows:

1\. Including Pluto’s haze emission is consistent with all current
observations. The best-fit model from our retrieval with haze yields about
1.63 mJy haze flux at 24 $\mu$m. In this case, Pluto’s haze emission could be
significant in the mid-infrared (e.g., 18–25.5 $\mu$m) compared with the
surface emission of Pluto or Charon.

2\. However, the haze flux has a large uncertainty due to the degeneracy of 24
$\mu$m thermal emission between Pluto’s haze and Charon’s surface. Even
without haze emission, the rotational lightcurves can be explained by the
surface emission from Pluto and a low-thermal-inertia Charon. Therefore, the
current data cannot constrain the haze flux well, ranging from 0 to 1.63 mJy,
with corresponding values of Charon’s diurnal thermal inertia from 5 to 23
MKS.

3\. The lightcurve shape can provide useful constraints on Pluto’s surface
properties. Our best-fit model indicates a diurnal thermal inertia of 8–10 MKS
on Pluto’s surface and a low emissivity of CH4 ice consistent with the small
grain case.

4\. JWST mid-infrared observations can separate Pluto from Charon and break
the degeneracy of emission between Pluto’s haze and Charon’s surface. We
predict that the rotational lightcurve of Pluto at 18 $\mu$m observed by JWST
in 2022-2023 will be relatively flat, which could directly assess the
existence of haze thermal emission as Pluto’s surface emission is only around
0.1 mJy assuming the same spectral emissivity at 18 $\mu$m as at 24 $\mu$m.

5\. The lightcurve variations of both Pluto and Charon in the JWST era range
from tens of $\mu$Jy at 18 $\mu$m to hundreds of $\mu$Jy at 25.5 $\mu$m. Their
variations can be characterized with JWST to constrain the thermal inertia and
ice albedo distributions on both bodies.

Previous atmospheric models that explained the cold temperature of Pluto
predict a wide range of haze emission at 24 $\mu$m. The haze flux is
approximately 2.48 mJy if its particles are composed of Titan-like tholins.
However, if the particles are composed of brighter ice particles, the flux
could decrease by more than an order of magnitude outside the ice absorption
band. The current observations cannot distinguish between the two scenarios,
and thus cannot constrain the cooling rate in Pluto’s atmosphere. At 18
$\mu$m, the haze flux is 1.5 (1.2) mJy for Titan-like tholins in 2004 (2023).
Outside the ice absorption band, the flux could be smaller than 0.1 mJy for
brighter ice particles. It is thus possible that JWST observations could help
to distinguish between the two scenarios and provide more information about
the cooling rate in Pluto’s atmosphere. Specifically, the mid-IR spectra from
JWST should be useful in identifying the solid-state features of ice
particles, thereby aiding the differentiation between the two potential
scenarios.

It is important to note the assumptions made in this study. First, we assume
that Pluto’s surface map has not changed significantly in the past two decades
when calculating the lightcurves of Spitzer in 2004 and Herschel in 2012
(Figure 3), and JWST in 2023 (Figure 5), based on the derived map from New
Horizons in 2015 (Figure 1). Although a time-invariant map in our model can
explain all observations, Pluto’s surface is expected to change over time. As
the surface maps of Pluto and Charon are crucial inputs when modeling
rotational emission lightcurves, a future study should include a seasonal
change of the surface ice and albedo distribution. That would require a
combined analysis of the New Horizons data revealing the limb-darkening
information, multi-decadal ground-based observations of Pluto and Charon
(Buratti et al., 2003, 2015), and an atmosphere-surface evolution model (Toigo
et al., 2015; Forget et al., 2017; Young, 2017; Bertrand et al., 2020).

Second, emissivity and thermal inertia are also important factors in
determining surface temperature besides albedo. Due to limited information, we
follow the approach of Lellouch et al. (2011) and divide Pluto’s surface into
three component units based on the albedo map and assume the same emissivity
within each unit. However, recent principal component analysis using the LEISA
observations has revealed several more ice units (Schmitt et al., 2017;
Gabasova et al., 2021; Emran et al., 2023). Different mixtures of ice could
have varying porosity and grain size, which would affect emissivity (Labed &
Stoll, 1991; Stansberry et al., 1996a). Additionally, we assume the same
thermal inertia for the entire globe, but it could vary with location and
depth due to changes in composition. Future work with better data may relax
these assumptions to retrieve the thermal inertia and emissivity of individual
geological units and connect them with the ice compositions from spectroscopic
measurements.

Finally, our model only accounts for the diurnal variation of surface and
subsurface temperatures, as we focus on rotational lightcurves within a day.
We assume that the heat flux from the deepest substrate in our model is zero,
which is not entirely accurate due to the large thermal inertia in the deep
layers of Pluto on a seasonal timescale. Although this assumption seems
working well for the infrared data, seasonal changes have been observed in
radio observations to probe the deep subsurface (e.g., Gurwell et al., 2011;
Butler et al., 2017; Bird et al., 2019). Future work should consider both
diurnal and seasonal changes like Lellouch et al. (2016), and combine both
infrared and radio data to finally understand the short-term and long-term
evolution of the surface, subsurface, and atmospheric haze on Pluto.

##  5 Acknowledgments

This work is supported by NASA Solar System Workings Grant 80NSSC19K0791. We
thank the reviewer for their valuable and constructive comments. We thank Dr.
Bonnie Buratti for the discussions about albedo data from New Horizons.
Retrievals are performed on lux supercomputer at UC Santa Cruz, funded by NSF
MRI grant AST 1828315.

## References

  * Bertrand et al. (2020) Bertrand, T., Forget, F., White, O., et al. 2020, Journal of Geophysical Research (Planets), 125, e06120. doi:10.1029/2019JE006120
  * Bertrand & Forget (2016) Bertrand, T. & Forget, F. 2016, Nature, 540, 86. doi:10.1038/nature19337
  * Bertrand & Forget (2017) Bertrand, T. & Forget, F. 2017, Icarus, 287, 72. doi:10.1016/j.icarus.2017.01.016
  * Bird et al. (2019) Bird, M. K., Linscott, I. R., Tyler, G. L., et al. 2019, Icarus, 322, 192. doi:10.1016/j.icarus.2019.01.004
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125. doi:10.1051/0004-6361/201322971
  * Buie & Grundy (2000) Buie, M. W. & Grundy, W. M. 2000, Icarus, 148, 324. doi:10.1006/icar.2000.6509
  * Buie et al. (1992) Buie, M. W., Tholen, D. J., & Horne, K. 1992, Icarus, 97, 211. doi:10.1016/0019-1035(92)90129-U
  * Buie et al. (1997) Buie, M. W., Tholen, D. J., & Wasserman, L. H. 1997, Icarus, 125, 233. doi:10.1006/icar.1996.5624
  * Buie et al. (2010a) Buie, M. W., Grundy, W. M., Young, E. F., et al. 2010, AJ, 139, 1117. doi:10.1088/0004-6256/139/3/1117
  * Buie et al. (2010b) Buie, M. W., Grundy, W. M., Young, E. F., et al. 2010, AJ, 139, 1128. doi:10.1088/0004-6256/139/3/1128
  * Buratti & Veverka (1983) Buratti, B. & Veverka, J. 1983, Icarus, 55, 93. doi:10.1016/0019-1035(83)90053-2
  * Buratti et al. (2017) Buratti, B. J., Hofgartner, J. D., Hicks, M. D., et al. 2017, Icarus, 287, 207. doi:10.1016/j.icarus.2016.11.012
  * Buratti et al. (2003) Buratti, B. J., Hillier, J. K., Heinze, A., et al. 2003, Icarus, 162, 171. doi:10.1016/S0019-1035(02)00068-4
  * Buratti et al. (2015) Buratti, B. J., Hicks, M. D., Dalba, P. A., et al. 2015, ApJ, 804, L6. doi:10.1088/2041-8205/804/1/L6
  * Buratti et al. (2019) Buratti, B. J., Hicks, M. D., Hillier, J. H., et al. 2019, ApJ, 874, L3. doi:10.3847/2041-8213/ab0bff
  * Butler et al. (2017) Butler, B. J., Gurwell, M. A., Lellouch, E., et al. 2017,
  * Cheng et al. (2008) Cheng, A. F., Weaver, H. A., Conard, S. J., et al. 2008, Space Sci. Rev., 140, 189. doi:10.1007/s11214-007-9271-6
  * Cheng et al. (2017) Cheng, A. F., Summers, M. E., Gladstone, G. R., et al. 2017, Icarus, 290, 112. doi:10.1016/j.icarus.2017.02.024
  * Emran et al. (2023) Emran, A., Dalle Ore, C. M., Ahrens, C. J., et al. 2023, , 4, 15. doi:10.3847/PSJ/acb0cc
  * Fan et al. (2022) Fan, S., Gao, P., Zhang, X., et al. 2022, Nature Communications, 13, 240. doi:10.1038/s41467-021-27811-6
  * Fan et al. (2019) Fan, S., Shemansky, D. E., Li, C., et al. 2019, Earth and Space Science, 6, 1057. doi:10.1029/2018EA000477
  * Forget et al. (2017) Forget, F., Bertrand, T., Vangvichith, M., et al. 2017, Icarus, 287, 54. doi:10.1016/j.icarus.2016.11.038
  * Fray & Schmitt (2009) Fray, N. & Schmitt, B. 2009, Planet. Space Sci., 57, 2053. doi:10.1016/j.pss.2009.09.011
  * Gabasova et al. (2021) Gabasova, L. R., Schmitt, B., Grundy, W., et al. 2021, Icarus, 356, 113833. doi:10.1016/j.icarus.2020.113833
  * Gladstone et al. (2016) Gladstone, G. R., Stern, S. A., Ennico, K., et al. 2016, Science, 351, aad8866. doi:10.1126/science.aad8866
  * Glasse et al. (2015) Glasse, A., Rieke, G. H., Bauwens, E., et al. 2015, PASP, 127, 686. doi:10.1086/682259
  * Grundy et al. (1999) Grundy, W. M., Buie, M. W., Stansberry, J. A., et al. 1999, Icarus, 142, 536. doi:10.1006/icar.1999.6216
  * Grundy et al. (2016) Grundy, W. M., Binzel, R. P., Buratti, B. J., et al. 2016, Science, 351, aad9189. doi:10.1126/science.aad9189
  * Gurwell et al. (2011) Gurwell, M., Butler, B., & Moullet, A. 2011, EPSC-DPS Joint Meeting 2011, 2011, 271
  * Hayne et al. (2017) Hayne, P. O., Bandfield, J. L., Siegler, M. A., et al. 2017, Journal of Geophysical Research (Planets), 122, 2371. doi:10.1002/2017JE005387
  * Hinson et al. (2017) Hinson, D. P., Linscott, I. R., Young, L. A., et al. 2017, Icarus, 290, 96. doi:10.1016/j.icarus.2017.02.031
  * Hofgartner et al. (2023) Hofgartner, J. D., Buratti, B. J., Beyer, R. A., et al. 2023, , 4, 132. doi:10.3847/PSJ/ace3ab
  * Holler et al. (2017) Holler, B. J., Young, L. A., Buie, M. W., et al. 2017, Icarus, 284, 394. doi:10.1016/j.icarus.2016.12.003
  * Jovanović et al. (2021) Jovanović, L., Gautier, T., Broch, L., et al. 2021, Icarus, 362, 114398. doi:10.1016/j.icarus.2021.114398
  * Jovanović et al. (2020) Jovanović, L., Gautier, T., Vuitton, V., et al. 2020, Icarus, 346, 113774. doi:10.1016/j.icarus.2020.113774
  * Khare et al. (1984) Khare, B. N., Sagan, C., Arakawa, E. T., et al. 1984, Icarus, 60, 127. doi:10.1016/0019-1035(84)90142-8
  * Labed & Stoll (1991) Labed, J. & Stoll, M. P. 1991, International Journal of Remote Sensing, 12, 2299. doi:10.1080/01431169108955259
  * Lavvas et al. (2021) Lavvas, P., Lellouch, E., Strobel, D. F., et al. 2021, Nature Astronomy, 5, 289. doi:10.1038/s41550-020-01270-3
  * Lellouch et al. (2000) Lellouch, E., Laureijs, R., Schmitt, B., et al. 2000, Icarus, 147, 220. doi:10.1006/icar.2000.6441
  * Lellouch et al. (2016) Lellouch, E., Santos-Sanz, P., Fornasier, S., et al. 2016, A&A, 588, A2. doi:10.1051/0004-6361/201527675
  * Lellouch et al. (2013) Lellouch, E., Santos-Sanz, P., Lacerda, P., et al. 2013, A&A, 557, A60. doi:10.1051/0004-6361/201322047
  * Lellouch et al. (2011) Lellouch, E., Stansberry, J., Emery, J., et al. 2011, Icarus, 214, 701. doi:10.1016/j.icarus.2011.05.035
  * Line et al. (2011) Line, M. R., Vasisht, G., Chen, P., et al. 2011, ApJ, 738, 32. doi:10.1088/0004-637X/738/1/32
  * Line et al. (2013) Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137. doi:10.1088/0004-637X/775/2/137
  * Mackay (2003) Mackay, D. J. C. 2003, Information Theory, Inference and Learning Algorithms, by David J. C. MacKay, pp. 640. ISBN 0521642981. Cambridge, UK: Cambridge University Press, October 2003., 640
  * Meza et al. (2019) Meza, E., Sicardy, B., Assafin, M., et al. 2019, A&A, 625, A42. doi:10.1051/0004-6361/201834281
  * Moore et al. (2016) Moore, J. M., McKinnon, W. B., Spencer, J. R., et al. 2016, Science, 351, 1284. doi:10.1126/science.aad7055
  * Moran et al. (2022) Moran, S. E., Hörst, S. M., He, C., et al. 2022, Journal of Geophysical Research (Planets), 127, e06984. doi:10.1029/2021JE006984
  * Nimmo et al. (2017) Nimmo, F., Umurhan, O., Lisse, C. M., et al. 2017, Icarus, 287, 12. doi:10.1016/j.icarus.2016.06.027
  * Ohno & Zhang (2019) Ohno, K. & Zhang, X. 2019, ApJ, 874, 1. doi:10.3847/1538-4357/ab06cc
  * Protopapa et al. (2021) Protopapa, S., Cook, J. C., Grundy, W. M., et al. 2021, The Pluto System After New Horizons, 433. doi:10.2458/azu_uapress_9780816540945-ch019
  * Reuter et al. (2008) Reuter, D. C., Stern, S. A., Scherrer, J., et al. 2008, Space Sci. Rev., 140, 129. doi:10.1007/s11214-008-9375-7
  * Schmitt et al. (2017) Schmitt, B., Philippe, S., Grundy, W. M., et al. 2017, Icarus, 287, 229. doi:10.1016/j.icarus.2016.12.025
  * Skilling (2004) Skilling, J. 2004, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, 735, 395. doi:10.1063/1.1835238
  * Spencer et al. (1989) Spencer, J. R., Lebofsky, L. A., & Sykes, M. V. 1989, Icarus, 78, 337. doi:10.1016/0019-1035(89)90182-6
  * Stansberry et al. (1996a) Stansberry, J. A., Pisano, D. J., & Yelle, R. V. 1996, Planet. Space Sci., 44, 945. doi:10.1016/0032-0633(96)00001-3
  * Stansberry et al. (1996b) Stansberry, J. A., Spencer, J. R., Schmitt, B., et al. 1996, Planet. Space Sci., 44, 1051. doi:10.1016/0032-0633(96)00084-0
  * Stern et al. (2017) Stern, S. A., Kammer, J. A., Gladstone, G. R., et al. 2017, Icarus, 287, 124. doi:10.1016/j.icarus.2016.09.019
  * Stern et al. (2021) Stern, S. A., White, O. L., McGovern, P. J., et al. 2021, Icarus, 356, 113805. doi:10.1016/j.icarus.2020.113805
  * Stern et al. (1997) Stern, S. A., Buie, M. W., & Trafton, L. M. 1997, AJ, 113, 827. doi:10.1086/118304
  * Strobel et al. (1996) Strobel, D. F., Zhu, X., Summers, M. E., et al. 1996, Icarus, 120, 266. doi:10.1006/icar.1996.0050
  * Strobel & Zhu (2017) Strobel, D. F. & Zhu, X. 2017, Icarus, 291, 55. doi:10.1016/j.icarus.2017.03.013
  * Toigo et al. (2015) Toigo, A. D., French, R. G., Gierasch, P. J., et al. 2015, Icarus, 254, 306. doi:10.1016/j.icarus.2015.03.034
  * Wan et al. (2021) Wan, L., Zhang, X., & Bertrand, T. 2021, ApJ, 922, 244. doi:10.3847/1538-4357/ac25f2
  * Yelle & Lunine (1989) Yelle, R. V. & Lunine, J. I. 1989, Nature, 339, 288. doi:10.1038/339288a0
  * Young & Binzel (1993) Young, E. F. & Binzel, R. P. 1993, Icarus, 102, 134. doi:10.1006/icar.1993.1038
  * Young et al. (1999) Young, E. F., Galdamez, K., Buie, M. W., et al. 1999, AJ, 117, 1063. doi:10.1086/300722
  * Young et al. (2001) Young, E. F., Binzel, R. P., & Crane, K. 2001, AJ, 121, 552. doi:10.1086/318008
  * Young (2017) Young, L. A. 2017, Icarus, 284, 443. doi:10.1016/j.icarus.2016.07.021
  * Young et al. (2018) Young, L. A., Kammer, J. A., Steffl, A. J., et al. 2018, Icarus, 300, 174. doi:10.1016/j.icarus.2017.09.006
  * Zhang et al. (2013) Zhang, X., West, R. A., Banfield, D., et al. 2013, Icarus, 226, 159. doi:10.1016/j.icarus.2013.05.020
  * Zhang et al. (2017) Zhang, X., Strobel, D. F., & Imanaka, H. 2017, Nature, 551, 352. doi:10.1038/nature24465

##  Appendix A Validation of the self-rotating thermophysical model against an
orbit-rotating model

We compare our self-rotating thermophysical model used in the retrieval to
that with the orbital revolution. Our simplified model includes only the
shallow subsurface, the planetary self-rotation, and assuming zero flux from
the deep subsurface. The orbit-rotating model, on the other hand, includes
both orbital revolution and self-rotation of the planetary body. Besides the
shallow subsurface, there is another deep subsurface where diurnal changes of
incoming solar flux could hardly have an impact, as described in Section 2.2.
Even though such simplification was suggested by Spencer et al. (1989), a
detailed comparison to the orbit-rotating model has not been reported.

In this validation, we chose an idealized setup. The input parameters of the
orbit-rotating model with two layers are the same as described in Section 2.2:
$\Gamma_{day}$ = 20 MKS till the depth of 0.07 m and additional
$\Gamma_{year}$ = 800 MKS underneath till 383.5 m. Other parameters of ice
properties are the same as the “free retrieval” scenario. We assume an orbit
similar to Pluto’s, which has an obliquity of 119.6∘, an eccentricity of
0.2444, and a semi-major axis of 39.48 AU. We first run the two-layer orbit-
rotating model using a coarse time step of 1 Pluto-day (6.4 Earth-days) for
several orbits so that the deep layer temperatures reach the steady state when
the deep temperatures at the same orbital position are consistent in two
successive orbits. Running three orbital periods seems to be enough for such a
long orbital period of 14179 self-rotating days from our test. Then we select
a starting point at the orbit as the initial condition (e.g., perihelion in
1989), and run simulations with both orbital motion and self-rotation, in a
finer time step of 1/1800 Pluto-day. Here we choose 800, 1600, and 2400 Pluto-
days after perihelion as the final positions on Pluto’s orbit for calculating
lightcurves (blue in Figure A1). We also run our simplified one-layer model
only with self-rotation at those final positions. After ensuring the diurnal
temperature variation reaches a steady state, the self-rotating lightcurves
(orange) are compared to orbit-rotating ones. The difference between these two
models (green) is smaller than 5% in 2003 and less than 1% in 2017 and 2031.
As our retrieval requires a fast forward model, this difference is acceptable,
trading off the accuracy and the efficiency.

Figure A1: Comparison of Pluto’s rotational lightcurves from the orbit-
rotating model (blue) and self-rotating model (orange) at 24 $\mu$m. Their
difference (green) is less than 5% of the total emission flux in 2003 and
smaller than 1% in 2017 and 2031.