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

# Effects of Haze Radiation and Eddy Heat Transport on the Thermal Structure of Pluto's Lower Atmosphere

Linfeng Wan (万霖丰)<fnref n="1" />, Xi Zhang<fnref n="1" />, and Tanguy Bertrand<fnref n="2" />,<fnref n="3" />

Received 2021 July 1; revised 2021 September 10; accepted 2021 September 10; published 2021 December 3

## Abstract

The temperature profile of Pluto's atmosphere has generally been assumed in a radiative-conductive equilibrium. Recent studies further highlighted the importance of radiative heating and cooling effects by haze particles. In this study, we update results from Zhang et al. by taking into account the icy haze composition proposed by Lavvas et al., and find that radiation of such an icy haze could still dominate the energy balance in the middle and upper atmosphere and explain the cold temperature observed by New Horizons. However, additional considerations are needed to explain the rapid decrease in temperature toward the icy surface at altitudes <25 km. We propose that vertical eddy heat transport might help maintain radiative-diffusive equilibrium in the lower atmosphere. In this scenario, our radiative-conductive-diffusive model (including both gas and haze) would match observations if the eddy diffusivity is on the order of $10^{-3} \text{ cm}^2 \text{ s}^{-1}$. Alternatively, if eddy heat transport is not effective on Pluto, in order to match observations, haze albedo must increase rapidly with decreasing altitude and approach unity near the surface. This is a plausible result of additional ice condensation and/or cloud formation. In this scenario, haze radiation might still dominate over gas radiation and heat conduction to maintain radiative equilibrium. Better constraints on haze albedo at ultraviolet and visible wavelengths would be a key to distinguish these two scenarios. Future mid-infrared observations from the James Webb Space Telescope could also constrain the thermal emission and haze properties in Pluto's lower atmosphere.

Unified Astronomy Thesaurus concepts: Pluto (1267)

## 1. Introduction

The New Horizons spacecraft flew by Pluto in 2015, providing a unique opportunity to closely observe its intriguing atmosphere and icy surface (Stern et al. 2015). Globally present haze is clearly seen with separate layers extending up to hundreds of kilometers. Pluto's atmospheric temperature and composition have also been derived from the New Horizons observations (Gladstone et al. 2016; Cheng et al. 2017; Young et al. 2018). From the radio occultation measurements, Pluto's atmospheric temperature over the Sputnik Planitia (during the ingress phase) ramps up rapidly from ~38.9 K near the surface to ~107 K at 25 km (Hinson et al. 2017). Ultraviolet occultation observations further constrained the temperature in the middle and upper atmosphere (above 25 km). It was found that the temperature decreases gradually from 25 to 475 km in altitude and remains almost isothermal of ~70 K above 475 km (Young et al. 2018; Gladstone & Young 2019). The temperature profile derived from the egress data exhibits a similar vertical structure (Hinson et al. 2017; Young et al. 2018).

The observed temperature in the middle and upper atmosphere is over 30 K lower than the theoretical prediction before the New Horizons encounter (Zhu et al. 2014). Strobel & Zhu (2017) proposed that water vapor could be a significant coolant, which is not likely to be possible, since the required amount would be supersaturated by many orders of magnitude. Zhang et al. (2017) argued that organic haze, rather than gases like H₂O or HCN, is a more reasonable coolant to explain Pluto's cold temperature in the middle and upper atmosphere. The thermal cooling and solar heating rates of haze are approximately equal, both of which are about 100 times greater than those of gases. Recently, Lavvas et al. (2021) used a sophisticated microphysical model to investigate the formation and composition of Pluto's haze. They proposed that

Pluto's haze may be primarily made of HCN ice above ~350 km and hydrocarbon ices such as C₄H₂ below 350 km. This is very different from the haze composition assumed in Zhang et al. (2017), which only considered HCN ice and laboratory Titan haze analogs instead of hydrocarbon ices. Using optical properties of the C₄H₂ ice, Lavvas et al. (2021) estimated the thermal cooling rate at 200 km is ~20 times smaller than that of Titan haze analogs in Zhang et al. (2017). On the other hand, icy haze might have a larger single scattering albedo than Titan haze analogs at visible wavelengths, resulting in a smaller solar heating rate. To date, Pluto's haze albedo and optical constants have not been constrained by the New Horizons data, and the radiative effects of icy haze on Pluto's temperature have not been quantitatively evaluated.

Including the radiative effect of haze might also challenge our current understanding of the thermal structure of Pluto's lower atmosphere below 25 km. Previously, it was believed that gas radiation and heat conduction are balanced in the lower atmosphere, leading to a positive temperature gradient with altitude (Gladstone & Young 2019). It has been suggested that haze radiation is at least one order of magnitude stronger than gas radiation and heat conduction (Zhang et al. 2017), so the haze radiative effect should also dominate the lower atmosphere. This implies that Pluto's lower atmosphere might be in a radiative equilibrium maintained by haze rather than the previously assumed radiative-conductive-equilibrium. However, there is a dilemma. Because the observed temperature decreases rapidly from its maximum (over 100 K) at 25 km to its minimum (less than 40 K) near the surface, the haze thermal cooling rate must decrease rapidly as well. But the haze solar heating rate should be almost constant if haze albedo stays roughly unchanged with height. As a result, the vertical profiles

<footnote n="1">¹Department of Earth and Planetary Sciences, University of California Santa Cruz, CA 95064, USA; wanmson@gmail.com</footnote>

<footnote n="2">²LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195, Meudon, France</footnote>

<footnote n="3">³NASA Ames Research Center, Space Science Division, Moffett Field, CA 94035, USA</footnote>
of haze heating and cooling cannot balance each other at each altitude in the lower atmosphere, where haze cooling becomes much weaker than haze heating near the surface. As we will show later, this heat imbalance will lead to overestimated temperature relative to the observed temperature under a radiative equilibrium.

To resolve this problem, in this study we propose and examine two possible solutions. In the first scenario, atmospheric dynamical mixing can transport heat and act as a cooling effect in the lower atmosphere. This vertical eddy transport process also occurs in Triton's lower atmosphere (Yelle et al. 1991) and other planetary atmospheres (e.g., Etling & Brown 1993; Forget et al. 1999; Colaïtis 2013). Alternatively, in the second scenario, haze heating may decrease near the surface if haze albedo increases toward the surface, and haze heating always balances haze cooling in the lower atmosphere.

We develop a 1D radiative–conductive–diffusive model including the previously ignored eddy heat transport to investigate physical processes and quantify possible scenarios in Pluto’s lower atmosphere. We introduce our model in Section 2, followed by the simulated temperature profiles for different scenarios in Section 3. Finally, this paper is concluded in Section 4.

## 2. Radiative–Conductive–Diffusive Model

Based on a previously published radiative–conductive model that includes heat conduction, adiabatic cooling, and radiation of both gas and haze (Zhang et al. 2017), we consider an additional energy transfer process by eddy heat transport. The heat transfer equations for gas (Equation (1a)) and haze (Equation (1b)) are:

$$ \rho_g c_{p,g} \frac{\partial T_g}{\partial t} = R_{\text{rad,g}} + R_{\text{cond}} + R_{\text{adb}} + R_{\text{eddy}} + R_{\text{pg}}, \quad (1a) $$

$$ \rho_p c_{p,p} \frac{\partial T_p}{\partial t} = R_{\text{rad,p}} - R_{\text{pg}}, \quad (1b) $$

where gas molecules and haze particles are denoted as g and p in subscripts, respectively. $\rho$ (kg m$^{-3}$) represents the density, $c_p$ (J K$^{-1}$ kg$^{-1}$) is the specific heat, and T (K) is the temperature. All of the energy transfer terms R are in units of W m$^{-3}$. $R_{\text{rad}}$ is the net radiative effect of solar heating and thermal cooling. $R_{\text{cond}}$, $R_{\text{adb}}$, and $R_{\text{eddy}}$ are the gas heat conduction, the adiabatic cooling due to gas escape, and the eddy heat diffusion, respectively. $R_{\text{pg}}$ is the heat exchange between gas molecules and haze particles due to collisions. The latent heating due to hydrocarbon/nitrile ice condensation is much smaller than the solar heating and is neglected in our model.

Pluto’s thin atmosphere is primarily in the gas phase, such that $\rho_g c_{p,g} \gg \rho_p c_{p,p}$. Also, gas molecules and haze particles are in the thermal equilibrium (i.e., $T_g=T_p$) due to efficient collisions (Zhang et al. 2017). Under these assumptions, we can combine Equations (1a) and (1b) into a single equation of gas heat transfer:

$$ \rho_g c_{p,g} \frac{\partial T_g}{\partial t} = R_{\text{rad,g}} + R_{\text{rad,p}} + R_{\text{cond}} + R_{\text{adb}} + R_{\text{eddy}}. \quad (2) $$

The radiative heating and cooling terms of gas ($R_{\text{rad,g}}$) and haze ($R_{\text{rad,p}}$) will be discussed in Sections 2.1 and 2.2, respectively. The heat conduction term ($R_{\text{cond}}$) and the adiabatic cooling term ($R_{\text{adb}}$) are the same as those in Strobel et al. (1996) and Zhang et al. (2017), whose expressions are shown in Appendix. The eddy heat

<figure><img src="image_2.png" /><figcaption>Figure 1. Volume mixing ratio profiles of Pluto’s atmospheric gases, assumed for our model and based on observations from New Horizons (N2, CH4, and C2H2; Gladstone et al. 2016; Young et al. 2018) and ALMA (CO and HCN; Lellouch et al. 2017). CO, HCN, and C2H2 are increased by factors of 2 × 10$^3$, 10$^4$, and 10$^2$ for illustrative purposes, respectively.</figcaption></figure>

diffusion term, $R_{\text{eddy}}$, can be represented as:

$$ R_{\text{eddy}} = \frac{1}{r^2} p^{\kappa} \frac{\partial}{\partial r} \left( r^2 \rho_g c_{p,g} K_{zz} p^{-\kappa} \left( \frac{\partial T_g}{\partial r} + \frac{g}{c_{p,g}} \right) \right), \quad (3) $$

where r (m) is the distance to Pluto’s center, p (Pa) is the pressure, g (m s$^{-2}$) is the gravity, and $\kappa$ is defined as the ratio of the universal gas constant and the specific heat. $K_{zz}$ (m$^2$ s$^{-1}$) is the eddy heat diffusivity, a key parameter in our model to explain the lower atmospheric temperature.

We solve Equation (2) by using the Crank-Nicolson method. The atmosphere below ~1800 km is divided into 176 uneven layers, whose thickness increases with height by a factor of 1.03 between adjacent layers (Strobel et al. 1996). In this study, we consider the temperature profile over the Sputnik Planitia as an example. Thus the lower boundary condition of temperature is fixed as 38.9 K. The upper boundary condition is assumed isothermal with no heat flux. More numerical details of the model are in Appendix.

### 2.1. Gas

Gases in Pluto’s atmosphere are mainly molecular nitrogen (N2) and methane (CH4), with minor constituents of acetylene (C2H2), carbon monoxide (CO), and hydrogen cyanide (HCN). Vertical profiles of gas abundances in Figure 1 are obtained from observations by New Horizons (N2, CH4, and C2H2; Gladstone et al. 2016; Young et al. 2018) and ALMA (CO and HCN; Lellouch et al. 2017). The mole fraction of CO, 5 × 10$^{-4}$, is assumed constant with altitude. The CH4 mixing ratio is assumed as 0.3% near the surface (Young et al. 2018) and is well mixed in the lower atmosphere where it is not strongly affected by photochemical destruction (Strobel & Zhu 2017). But the CH4 mole fraction increases with altitude due to its molecular separation from N2 above the homopause where molecular diffusion of the species exceeds eddy diffusion.

The gas radiative heating and cooling rates are calculated by using the Correlated-k method in Strobel et al. (1996) and Zhang et al. (2017). We consider the solar heating of CH4 in the infrared vibrational bands at 1.7, 2.3, 3.3, and 7.6 µm under
the non-local thermodynamic equilibrium (non-LTE), as well as the heating of CH₄ and N₂ in the far-ultraviolet and extreme-ultraviolet. For the gas thermal cooling, we consider rotational bands of CO and HCN under the LTE, the vibrational band of CH₄ at 7.6 µm, and that of C₂H₂ at 13.7 µm under the non-LTE.

## 2.2. Haze

Pluto's haze is clearly seen in the New Horizons images below 350 km, and its extinction coefficient is measured from the solar ultraviolet occultation (Gladstone et al. 2016; Cheng et al. 2017; Young et al. 2018). Based on previous studies of Titan's hydrocarbon haze (e.g., Lavvas et al. 2013), photo-chemistry and ionchemistry of N₂, CO, and CH₄ might also occur in Pluto's upper atmosphere where nanometer-sized monomers could form and then settle down. These monomers will coagulate into micrometer-sized fractal aggregates (Bertrand & Forget 2017; Gao et al. 2017; Zhang et al. 2017). The latest microphysical modeling results suggested that Pluto's haze is formed differently from Titan's, which is likely to be ice forming through direct condensation of photochemical products (Lavvas et al. 2021). In this "icy haze" situation, HCN ice particles initially form in the upper atmosphere, settle down due to gravity, and grow through subsequent condensation by hydrocarbon ices and heterogeneous chemistry. The major composition of the simulated haze thus shifts from HCN ice above 350 km to C₄H₂ ice below 350 km.

Multiple ice mixtures will probably lead to different radiative energy fluxes from Titan haze analogs assumed in Zhang et al. (2017). The Titan-like organic haze primarily absorbs in the ultraviolet and visible and emits strongly in the infrared, leading to the large rates of radiative heating and cooling. Due to a large uncertainty of haze refractive indices, its heating and cooling rates can vary by one to two orders of magnitude but are still consistently larger than gas rates (Zhang et al. 2017). However, if Pluto's haze is mostly made of ice mixtures, the haze heating and cooling rates are likely to be smaller than those for the Titan-like organic haze. The ice heating rate could be smaller because of the higher albedo in the visible. The infrared cooling mechanism by ices is fundamentally different from that by tholin-like particles because the ices can only emit at their specific vibrational bands in the infrared while the tholin emission can occur across a broad continuum. Lavvas et al. (2021) estimated the infrared cooling rate of the C₄H₂ ice, and their "preliminary evaluation of the cooling rate based only on the narrow absorption bands of C₄H₂ ice near 40 µm suggests a value of 10⁻¹⁰ W m⁻³ at 200 km that is ~20 times lower than the previous estimate based on Titan's analogues." If this were true, is the icy haze still able to explain the observed temperature profile from New Horizons, despite smaller heating and cooling rates? Here we revisit the radiative energy balance in Pluto's atmosphere with the proposed icy haze.

To calculate icy haze radiation in our model, we first compare the haze properties in Lavvas et al. (2021) with those in Zhang et al. (2017), which are mainly based on Gao et al. (2017). A rough comparison shows that vertical distributions of the particle size in the two models are not very different. For example, the particle size is about 0.5 µm near the ground and 0.1 µm at around 100 km in both models. The vertical profiles of haze abundance are not that different either, as both are constrained by the Alice data from New Horizons. The main difference between icy haze and organic haze is the refractive indices, which are highly uncertain for both organic haze and

multiple ice mixtures. In particular, the haze heating and cooling rates critically depend on the imaginary part of the particle refractive index ($k$), which is very poorly known for hydrocarbon and nitrile ices (e.g., Khanna et al. 1988; Moore et al. 2010; Hudson et al. 2014). In the UV, a preliminary estimation of $k$ for C₄H₂ ice particles is about 10⁻² (Lavvas et al. 2021), which is comparable to the values of some laboratory tholins (e.g., Ramirez et al. 2002). However, $k$ might steeply decrease toward visible wavelengths (e.g., for CH₄, Martonchik & Orton 1994), so it is expected that the ice is brighter than Titan's tholins in the visible. In the near-infrared, both amorphous and crystalline HCN ices at low temperatures have a $k$ value of about 5 × 10⁻⁴ at 2 µm (Moore et al. 2010), which is comparable to that of Titan haze analogs in Khare et al. (1984). However, ices might become absorptive in their absorption bands (e.g., Hudson et al. 2014).

The infrared cooling mainly occurs in the thermal and far-infrared wavelengths. $k$ in this range has been measured for hydrocarbon ices such as C₄H₂ (Khanna et al. 1988), C₂H₂, C₂H₄, and C₂H₆ (Hudson et al. 2014), and nitrile ices such as HCN and HC₃N (Moore et al. 2010). Because Pluto's atmospheric temperature varies from 40 to 100 K, only ices with absorption bands between 30–100 µm can cause effective infrared cooling. For all ices simulated in Lavvas et al. (2021), the two most abundant components are C₄H₂ and HCN. The absorption band of the C₄H₂ ice at 40 µm and the HCN band at around 60 µm should be the most important cooling bands in Pluto's atmosphere. Because the HCN absorption band is much stronger than the C₄H₂ ice, previous calculations show that the cooling rate by pure HCN ice particles can be as large as that by Titan's tholins (Zhang et al. 2017), while 20 times larger than that by the pure C₄H₂ ice (Lavvas et al. 2021).

Moreover, the detailed structure of ice mixtures should also play an important role and add another level of complexity. Anderson et al. (2018) studied the organic ice in Titan's atmosphere and pointed out that spectra from ice mixtures could behave very differently from their pure components. Their Figure 12 demonstrates that the ice mixture from co-condensation of HCN and HC₃N shows different spectra from the layered ice. Co-condensation should be expected in Pluto's ice particle formation through the heterogeneous chemistry as well (see Figure 3 in Lavvas et al. 2021), but the resultant ice structure is unknown. Therefore, it is difficult to estimate the single scattering albedo and thermal infrared opacity for Pluto's ice particles proposed in Lavvas et al. (2021) based on our current knowledge of pure ices. In this study, we do not attempt a very detailed investigation. Instead, we solve the icy haze radiation using the multiple-scattering radiative transfer code in Zhang et al. (2017) as for the case of organic haze but with the following modifications to take into account the impact of the ice component of haze:

(1) Below 200 km, C₄H₂ dominates the ice composition (Lavvas et al. 2021). Because the pure C₄H₂ ice behaves less effectively than HCN at absorption bands, we just assume a pure C₄H₂ ice in the cooling calculation and assume it is the lower limit of the cooling rate by possible ice mixtures in Pluto's lower atmosphere. The reasoning behind this choice is that if this possibly underestimated haze cooling still dominates the energy balance in the atmosphere, the proposed explanation of Pluto's temperature profile using haze particles by Zhang et al. (2017) remains valid. In our calculation, we apply a factor of 1/20 on the infrared optical
<table><caption>**Table 1**
Eddy Diffusivity ($K_{zz}$) Derived from Several Chemical Models with the Near-Surface CH$_4$ Mixing Ratio ($f_{CH_4}$)</caption><thead><tr><th>$f_{CH_4}$</th><th>$K_{zz}$ (cm<sup>2</sup> s<sup>-1</sup>)</th><th>Homopause (km)</th><th>References</th></tr></thead><tbody><tr><td>0.60%-0.84%</td><td>7.5 × 10<sup>5</sup> at the surface; 3 × 10<sup>6</sup> at 210 km</td><td>390</td><td>Gladstone et al. (2016)</td></tr><tr><td>0.4%</td><td>1 × 10<sup>3</sup><br>1 × 10<sup>6</sup> under 150 km; 3 × 10<sup>6</sup> above 150 km</td><td></td><td>Wong et al. (2017)</td></tr><tr><td>0.28%-0.35%</td><td>5.5 × 10<sup>2</sup> - 4 × 10<sup>3</sup></td><td>2-12</td><td>Luspay-Kuti et al. (2017)</td></tr><tr><td>0.45%</td><td>3 × 10<sup>4</sup></td><td>110</td><td>Young et al. (2018)</td></tr><tr><td>0.4%</td><td>1 × 10<sup>3</sup> near the surface; 10<sup>3</sup>-10<sup>4</sup> above</td><td></td><td>Krasnopolsky (2020)</td></tr><tr><td></td><td></td><td></td><td>Lavvas et al. (2021)</td></tr></tbody></table>

depth calculated in Zhang et al. (2017) to mimic the C$_4$H$_2$ ice radiative effect estimated in Lavvas et al. (2021) and to obtain a thermal cooling rate of the haze that is reduced by about a factor of 20 compared to the organic haze case.

(2) Because the HCN ice still dominates the icy haze composition above 350 km (Lavvas et al. 2021), the HCN ice cooling rate in the upper atmosphere should still be comparable to that calculated in Zhang et al. (2017; see the HCN ice cooling rate in their Figure 2). We thus scale the infrared optical depth in Zhang et al. (2017) from 200 to 700 km by a factor of $(z/700)^{2.4}$, where z (km) is the altitude from surface. This scaling factor gradually varies from 1/20 at 200 km to unity at 700 km.

(3) Finally, we also apply the same scaling factor (1/20 below 200 km and increasing gradually to unity at 700 km) to the haze solar heating rate profile calculated in Zhang et al. (2017). This choice of scaling factor is to approximate the particle albedo increase due to the ice component and allows us to match the observed temperature profile of New Horizons. For reference, decreasing the haze solar heating rate by a factor of 20 roughly corresponds to increasing the haze albedo from 0.5 to 0.975 in the visible.

The simulated results and detailed discussion of the energy balance will be presented in Section 3.

## 2.3. Eddy Diffusivity

As mentioned earlier, eddy heat transport could be important in Pluto's lower atmosphere. It can be approximated as a diffusional process. Eddy diffusion is generally used to treat atmospheric dispersion and mixing of energy and chemical species macroscopically. Tracer transport and heat mixing occur due to katabatic winds, general circulations, wave breaking, and turbulence generated by the near-surface ice condensation and sublimation (Bertrand & Forget 2016; Forget et al. 2017; Bertrand et al. 2020a, 2020b; Jacobs et al. 2021). Zalucha et al. (2011) tested the convective mixing of potential temperature in a putative troposphere of Pluto to explain the observed occultation light curves, but New Horizons did not find a convective troposphere on Pluto. Instead, the existing temperature inversion with a positive gradient of 3 K km$^{-1}$ from 0 to 25 km suggests that Pluto's lower atmosphere is stably stratified without strong large-scale convection. Our hypothetical eddy diffusion, if exists, might be a result of atmospheric waves and large-scale circulation instead of convection.

In our 1D framework, eddy diffusivity ($K_{zz}$) is a key parameter for vertical eddy transport. Although here we mainly focus on $K_{zz}$ for the heat diffusion, which could be different from that for chemical tracer transport, the latter might provide some reference for the eddy heat diffusivity. One can empirically determine $K_{zz}$ (for tracer transport) by fitting the

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Eddy diffusivity ($K_{zz}$) derived from Pluto's GCM (blue) by calculating with potential temperature (solid) and CH$_4$ mixing ratio (dashed) in Equation (4). Molecular diffusion coefficient (purple dotted) is also plotted for comparison (Krasnopolsky & Cruikshank 1999). Other lines are $K_{zz}$ derived from chemical models listed in Table 1. The red lines with arrows indicate a range between $5.5 \times 10^2$ and $4 \times 10^3$ cm$^2$ s$^{-1}$ from Young et al. (2018). The orange line is approximately the upper limit from Krasnopolsky (2020) that still matches the homopause in the CH$_4$ observations.</figcaption></figure>

distribution of a chemical tracer such as CH$_4$ with a chemical-transport model. CH$_4$ is the second most abundant species in Pluto's atmosphere and is not greatly affected by photolysis loss in the near-surface due to self-shielding, and thus is approximately well mixed below the homopause. Based on the observed CH$_4$ vertical profile above 100 km from New Horizons, previous chemical-transport simulations have been used to constrain $K_{zz}$ on Pluto. Their results are summarized in Table 1 and plotted in Figure 2. Wong et al. (2017) first constrained $K_{zz} = 10^3$ cm$^2$ s$^{-1}$. Later on, Gladstone et al. (2016), Luspay-Kuti et al. (2017), Mandt et al. (2017), and Lellouch et al. (2017) obtained a much larger $K_{zz} \sim 10^6$ cm$^2$ s$^{-1}$. Young et al. (2018) updated the CH$_4$ profile using a joint analysis of the Alice CH$_4$ data and Alice and REX N$_2$ measurements from New Horizons. Their results supported the smaller $K_{zz} \sim 10^3$ cm$^2$ s$^{-1}$. Recently, Krasnopolsky (2020) found that $3 \times 10^4$ cm$^2$ s$^{-1}$ can still fit the CH$_4$ data from Young et al. (2018). The microphysical model in Lavvas et al. (2021) used $K_{zz} = 10^3$ cm$^2$ s$^{-1}$ near the surface and 1–10 times larger at higher altitudes. The disagreement among the previous studies primarily originated from a degeneracy between the assumed near-surface methane mixing ratio and the derived $K_{zz}$. The near-surface CH$_4$ mixing ratio ranges from 0.3% to 0.8% in these chemical models, resulting in a large uncertainty range of derived $K_{zz}$ from $10^3-10^6$ cm$^2$ s$^{-1}$ (Table 1).
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) Temperature profiles of Pluto's atmosphere from 1D model simulations with icy haze (red, in this study) and organic haze (purple; Zhang et al. 2017). Black lines are nominal (solid) and flux-constrained (dashed) analytic temperature profiles from New Horizons, and black error bars correspond to derived temperatures during ingress (diamond) and egress (circle) phases (Young et al. 2018). (b) An inset of panel (a) for the atmosphere below 120 km, with the New Horizons radio occultation (REX) measurements during ingress (diamond) and egress (circle) phases (Hinson et al. 2017). (c) Radiative heating (solid) and cooling (dashed) rates of gases (gray) and haze (green), and heat conduction (blue) based on the observed temperature profile and the assumed icy haze properties.</figcaption></figure>

In principle, $K_{zz}$ can also be derived from 3D simulations with tracer transport using chemical species or potential temperature. In this view, the 1D eddy transport is a global-average approximation of the 3D effect (Zhang & Showman 2018). One can define the eddy-mean decomposition by $A = \bar{A} + A'$, where A represents any quantity with its globally averaged quantity $\bar{A}$, and $A'$ called eddy is the deviation from the mean. To estimate the global mean vertical flux of a tracer, we take Equation (6) in Zhang & Showman (2018):

$$ \overline{w\chi'} \approx -K_{zz} \frac{\partial \bar{\chi}}{\partial z}, \quad (4) $$

where $w$ (m s⁻¹) is the vertical wind speed, $z$ (m) is the altitude, and $\chi$ here takes one specific quantity of a tracer from the model. Using the simulated distributions of potential temperature and CH₄ as tracers from Pluto's general circulation model (GCM; Bertrand & Forget 2016; Forget et al. 2017; Bertrand et al. 2020a, 2020b), one can derive $K_{zz}$ by Equation (4) based on the #1 simulation for 2015 in Bertrand et al. (2020b). $K_{zz}$ derived from GCM in Figure 2 by heat transport (blue solid) and tracer transport (blue dashed) are not necessarily the same. The derived $K_{zz}$ for heat transport is about $1.7 \times 10^3$ cm² s⁻¹ near the surface and increases to $10^4$ cm² s⁻¹ at higher altitudes. The derived $K_{zz}$ for CH₄ transport is generally larger, starting from $6 \times 10^3$ cm² s⁻¹ near the surface, exhibiting a non-monotonic shape until 25 km, and reaching about $10^6$ cm² s⁻¹ at higher altitudes.

Moreover, $K_{zz}$ is closely related to the atmospheric homopause, above which molecular diffusion will take control over eddy diffusion. From the latest CH₄ profile in Young et al. (2018), the CH₄ homopause should be below ~100 km, putting the upper limit of $K_{zz}$ for CH₄ transport at $3 \times 10^4$ cm² s⁻¹ (orange line in Figure 2). $K_{zz}$ should not be too large (e.g., $10^5$ cm² s⁻¹), otherwise the homopause will be too high. $K_{zz}$ for heat transport derived from Pluto's GCM (blue solid line in Figure 2) is also smaller than this upper limit. We will test different values of $K_{zz}$ in the temperature simulation (Section 3.1), from which the New Horizons observations will in turn constrain $K_{zz}$.

### 3. Two Possible Scenarios to Explain Pluto's Temperature Profile

We first show a “traditional” case without eddy diffusion, following previous studies in Zhu et al. (2014), Strobel & Zhu (2017), and Zhang et al. (2017). The purpose of this test is to examine whether the newly proposed icy haze from Lavvas et al. (2021) could still explain the cold upper atmosphere where the eddy diffusion is not important above the homopause. Using the estimated optical properties of icy haze in Section 2.2, we found that the icy haze model is still able to explain the temperature profile in Pluto's middle and upper atmosphere (Figure 3(a)). For comparison we also plotted the simulated temperature profile from the previous organic haze model in Zhang et al. (2017). Even though the thermal cooling rate of the icy haze is much smaller than that of the organic haze in Zhang et al. (2017), it is still much greater than that of gases. Note that we have also applied a reduction factor of the solar heating rate to approximate the albedo increase of the icy haze (see Section 2.2) so that the icy haze heating can balance its thermal cooling. Given the large uncertainties of albedo and optical properties of such an icy haze, the model in Figure 3(a) is by no means the only solution to explain Pluto's cold upper atmosphere. However, our test confirms that it is still possible to explain the New Horizons data with the recently proposed icy haze on Pluto (Lavvas et al. 2021). The icy haze test also qualitatively agrees with the conclusion in Zhang et al. (2017) that haze radiation controls the energy balance and haze infrared cooling helps maintain a cold temperature in Pluto's upper atmosphere.

However, we find that the simulated temperature with icy haze is much hotter than observations in the lower atmosphere (Figure 3b). This is also true for the organic haze case in Zhang et al. (2017), which mostly focused on explaining the cold upper atmosphere. To understand the reason, we calculate the energy fluxes by using the observed temperature and icy haze properties in Figure 3(c), and find that heating and cooling rates of icy haze cannot be balanced. This imbalance is caused by the rapidly decreasing temperature toward the icy surface so that the haze cooling rate drastically changes from $10^{-8}$ W m⁻³ at 25 km to $10^{-10}$ W m⁻³ near the surface. But the haze heating
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Temperature profiles of Pluto's atmosphere from 1D model simulations with icy haze and eddy heat transport (red). Black lines are nominal (solid) and flux-constrained (dashed) analytic temperature profiles from New Horizons, and black error bars correspond to derived temperatures during ingress (diamond) and egress (circle) phases (Young et al. 2018). (b) Radiative heating (solid) and cooling (dashed) rates of gases (gray) and haze (green), heat conduction (blue), and eddy heat transport (magenta) in the model. (c) An inset of panel (a) below 120 km, with the New Horizons REX measurements during ingress (diamond) and egress (circle) phases (Hinson et al. 2017). Simulations with different $K_{zz}$ are also plotted in the same colors of Figure 2. (d) An inset of panel (b) below 120 km.</figcaption></figure>

rate remains roughly constant with $10^{-8}$ W m$^{-3}$ under 25 km because the solar flux is not attenuated significantly in the optical-thin lower atmosphere. Integrated from the ground to 25 km, the total excess haze heat flux in the lower atmosphere is on the order of $10^{-4}$ W m$^{-2}$. Gas radiation and heat conduction are negligible compared with haze radiation, and thus cannot compensate for the haze net heating.

To explain the thermal structure of Pluto's lower atmosphere, we propose two possible solutions. The first one is called the “radiative-diffusive equilibrium” scenario, in which eddy heat transport acts as a cooling effect to remove the excess heating in Pluto's lower atmosphere. Because potential temperature increases with altitude in the lower atmosphere, eddy diffusion of potential temperature transports heat from a higher altitude to a lower altitude. In principle, the eddy transport could either heat or cool a local atmospheric layer, depending on the local heat flux divergence and the curvature of the potential temperature profile. For each layer in Pluto's lower atmosphere, we find that the eddy heat flux transported downward from the upper layers is smaller than that transported farther down to the lower layers, resulting in a local cooling effect that could balance the haze net heating. Alternatively, in the second scenario, if eddy heat transport is not effective in the

stably stratified lower atmosphere on Pluto, the lower atmosphere might still be in the radiative equilibrium maintained by haze radiative heating and cooling. However, haze albedo in this case must increase toward the surface to reduce the haze solar heating so that it could balance the haze thermal cooling under 25 km. We call this the “bright-haze” scenario. Note that the two above-mentioned scenarios are just end members and are not mutually exclusive. We investigate the two end-member scenarios separately in the following two subsections.

### 3.1. Radiative–Diffusive Equilibrium Scenario

For the eddy heat transport to be significant, the eddy heat flux convergence must be comparable to the haze net heating rate in the lower atmosphere. In our icy haze model, we find that a constant $K_{zz} = 4 \times 10^{-3}$ cm$^2$ s$^{-1}$ (red solid line in Figure 2, the upper limit from Young et al. 2018) is able to explain the New Horizons radio occultation (REX for short) measurements in Pluto's lower atmosphere (Figure 4(c), red). Eddy heat transport acts as a cooling process with a rate of $3 \times 10^{-8}$ W m$^{-3}$ at 5 km, one order of magnitude greater than haze radiative cooling, gas radiation, and heat conduction. It implies that eddy heat transport and haze radiative heating, rather than previously believed heat
conduction and gas radiation, dominate the thermal structure of Pluto's lower atmosphere. In other words, the atmosphere is in a radiative-diffusive equilibrium instead of a radiative-conductive equilibrium. In this case, the atmospheric eddies diffuse the excess haze radiative flux on the order of $10^{-4}$ W m$^{-2}$ down to the surface, which is negligible compared with the absorbed insolation on Pluto (on the order of $10^{-1}$ W m$^{-2}$) and therefore has no significant impact on icy surface temperature and/or on volatile sublimation/condensation surface rates. Above 25 km, because the atmospheric density drops quickly with increasing altitude, the efficiency of eddy heat transport decreases rapidly and becomes comparable to or even less than gas radiation. Therefore, eddy heat transport with $K_{zz} = 4 \times 10^3$ cm$^2$ s$^{-1}$ has a negligible contribution to the energy balance when compared with haze radiation in the middle and upper atmosphere (Figure 4b), where our simulated temperature profile is similar to that in the icy model without eddy diffusion (red line in Figure 3).

The temperature observations from New Horizons also provide an additional constraint on $K_{zz}$ in the lower atmosphere. We test different values of $K_{zz}$ in our model and find that eddy cooling from the case with $K_{zz} = 3 \times 10^4$ cm$^2$ s$^{-1}$ (orange line in Figure 2; Krasnopolsky 2020) is too strong to explain the observations (Figure 4(c), orange). We also test $K_{zz}$ derived from Pluto's GCM using potential temperature as a tracer (blue solid line in Figure 2), which is within the $K_{zz}$ range from Young et al. (2018) below 15 km and is a factor of three smaller than its upper limit below 5 km (red dashed and solid lines in Figure 2). In this case, eddy heat transport is not effective at completely removing the excess haze heat, and the simulated temperature is still 5 K hotter at 10 km compared with the occultation data (Figure 4(c), blue). We also test the $K_{zz}$ derived from GCM using the CH$_4$ as a tracer (blue dashed line in Figure 2). The simulated lower atmospheric temperature (not shown here) is close to that using $K_{zz} = 3 \times 10^4$ cm$^2$ s$^{-1}$ (Figure 4(c), orange) and thus is too cold to explain the data. The simulated results above suggest that an eddy heat diffusivity of $K_{zz}$ is on the order of $10^3$ cm$^2$ s$^{-1}$ in Pluto's lower atmosphere for the icy haze albedo and opacities assumed in our model. It should also be emphasized that derived $K_{zz}$ strongly depend on albedo and opacities of the icy haze, which are very uncertain at this moment. Nevertheless, in this radiative-diffusive scenario, eddy heat transport seems necessary to explain the thermal structure of Pluto's lower atmosphere.

## 3.2. Bright Haze Scenario

If $K_{zz}$ is as small as $550$ cm$^2$ s$^{-1}$ (red dashed line in Figure 2, the lower limit from Young et al. 2018), eddy heat transport does not effectively remove the excess haze heating, and the simulated temperature is hotter than the observations by about 10 K at 10 km in our icy haze model (Figure 4(c), red dashed). If this were the realistic eddy diffusivity on Pluto, here we propose another solution to explain the temperature profile if icy haze albedo significantly increases toward the surface. To elaborate on this “bright-haze” scenario, we consider an “end-member” case without any eddy diffusion. The bright icy particles might still be able to maintain the radiative equilibrium instead of the radiative-diffusive equilibrium in Pluto's lower atmosphere. A possible cause of this albedo change is the additional ice coating. As temperature drops rapidly toward the cold icy surface, more ice might condense onto haze particles and increase its albedo (Lavvas et al. 2021). Whether this could occur in the real Pluto atmosphere is debatable. Since Pluto's haze particles are made of

condensed ice as proposed in Lavvas et al. (2021), icy haze could be bright already in the upper atmosphere. It seems difficult to further increase haze albedo in the lower atmosphere due to additional ice condensation. Another possible cause of the albedo change is cloud formation near the surface. If the temperature is cold enough, the icy haze could serve as cloud condensation nuclei for supersaturated methane and nitrogen over the icy surface (Forget et al. 2017; Stern et al. 2017; Bertrand et al. 2020b); although, the bright methane and nitrogen clouds on Pluto would be sparse anyway, and they still lack definitive evidence so far.

Hypothetically, if additional ice coating by condensation or cloud formation in the lower atmosphere would increase haze albedo and decrease the haze solar heating rate, it is possible for haze heating to balance infrared cooling. In the near-surface layers, since the infrared cooling rate is very small (Figure 5), the haze heating rate should also approach zero, suggesting that the albedo of ice-coated particles must reach almost unity near the surface. This is reminiscent of very bright haze with unity albedo that was detected in the lower atmosphere of Neptune's moon Triton during the Voyager flyby (Hillier et al. 1991).

To test this scenario, we set $K_{zz} = 0$ in our icy haze model and apply another scaling factor on the haze heating rate profile in the lower atmosphere to mimic the albedo increase. The multiplication factor decreases gradually from 1.0 at ~20 km to 0.003 near the surface. Correspondingly, haze albedo in the visible increases approximately from 0.975 at ~20 km to very close to unity near the surface. With the modified solar haze heating rate, the observed temperature can be reproduced (Figure 5). In this “bright-haze” scenario, haze heating and cooling are still much stronger than gas radiation and heat conduction (except very close to the surface) so that the lower atmosphere is in a radiative equilibrium regulated by the haze radiation.

In an extreme limit, if haze is even brighter than the assumed values above, the haze heating and cooling can be even weaker than gas radiation and heat conduction, and the lower atmosphere will be controlled by gas net heating and heat conduction. This goes back to the previous gas-only studies (e.g., Strobel & Zhu 2017). However, this extreme case requires very bright icy haze with very low infrared opacity (or high infrared albedo). Future characterization of haze albedo and its vertical profile near the surface is a key to understanding the details of energy balance in Pluto's lower atmosphere.

## 4. Summary and Discussions

In this paper, we revisit the physical processes that influence the thermal structure of Pluto's lower atmosphere. The problem originates from the large haze radiative effect that was proposed to explain the temperature profile of Pluto's middle and upper atmosphere (Zhang et al. 2017). Our calculation shows that even if Pluto's haze is mainly composed of HCN and hydrocarbon ices, haze radiation is still stronger than that of gas. However, if haze radiation still dominates in the lower atmosphere, the absorbing haze leads to greater heating near the cold surface than infrared cooling and cannot explain the temperature observations.

We propose two scenarios to overcome this energy imbalance problem. In the “radiative-diffusive equilibrium” scenario, eddy heat transport provides an effective cooling mechanism to balance the excess heating of haze. The required eddy diffusivity is on the order of $10^3$ cm$^2$ s$^{-1}$, which is
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Temperature profiles of Pluto's atmosphere from 1D model simulations with brighter icy haze below 20 km but no eddy heat transport (red). Black lines are nominal (solid) and flux-constrained (dashed) analytic temperature profiles from New Horizons, and black error bars correspond to derived temperatures during ingress (diamond) and egress (circle) phases (Young et al. 2018). (b) Radiative heating (solid) and cooling (dashed) rates of gases (gray) and haze (green), heat conduction (blue), and eddy heat transport (magenta) in the model. (c) An inset of panel (a) below 120 km, with the New Horizons REX measurements during ingress (diamond) and egress (circle) phases (Hinson et al. 2017). (d) An inset of panel (b) below 120 km.</figcaption></figure>

comparable with the chemical tracer diffusivity from chemical and dynamical models. On the other hand, if eddy diffusion is not strong enough to transport heat in the stably stratified lower atmosphere, we suggest that bright haze near the surface might greatly reduce its solar heating rate and thus explain the cold near-surface temperature on Pluto. The key of this “bright-haze” scenario is the single scattering albedo increasing from ~20 km to the ground, which might be a result of ice coating or cloud formation in the lower atmosphere. The near-surface haze albedo could approach as large as unity in the visible, which is similar to what has been found on Triton.

The two suggested scenarios might not be mutually exclusive. One could have both eddy heat transport and bright haze affecting the energy balance in Pluto's lower atmosphere. However, the large uncertainty of haze optical properties prevents us from doing a rigorous and quantitative investigation of these hypotheses.

Future efforts can be put forward in several directions. Haze made of different chemical compositions could exhibit different infrared emission features, for example, broad continuum for Titan's tholins (Zhang et al. 2017) but sharp spectral features for hydrocarbon ices (Lavvas et al. 2021). Therefore, infrared

observations by JWST would be useful to constrain the infrared properties and the cooling rate of Pluto's haze (Zhang et al. 2017). However, JWST data will not be sufficient to fully distinguish the two proposed scenarios, in which the thermal infrared properties of icy haze are assumed to be the same in both eddy-diffusion and bright-haze calculations in this study, leading to similar infrared emission spectra. Instead, shortwave observations of Pluto's haze are the key to constrain haze albedo and the solar heating rate in the ultraviolet and visible. One possible way is to directly measure the extinction coefficients from occultation at visible wavelengths. In combination with the reflectivity observations from New Horizons, the single scattering albedo of haze particles might be constrained, although the current New Horizons data are limited within several channels and several phase angles. Thus, the albedo and the scattering phase function of haze particles might require future observations to test the possible haze albedo change under 25 km in the “bright-haze” scenario. Moreover, there is a lack of laboratory measurements of the refractive indices of Pluto icy haze analogs. Some recent efforts have been put forward to constrain the chemical composition and optical constants of ice-free Pluto aerosol analogs
(Jovanović et al. 2020, 2021). However, similar measurements for icy particles are still lacking. Future laboratory investigation of ice-coated haze in a Pluto-like environment could also shed light on the optical properties and the radiative effect of Pluto's haze in the lower atmosphere. Finally, it should be kept in mind that the observed temperature profile in the near-surface atmosphere could vary from place to place depending on local conditions, as evidenced by the differences between the ingress and egress profiles. If neither of these profiles represents a good average for Pluto—which needs to be verified by future data— our proposed mechanisms are only related to local processes rather than global phenomena.

We thank Leslie Young for providing the derived CH₄ profile from Young et al. (2018). This work is supported by NASA Solar System Workings grant 80NSSC19K0791.

## Appendix
### Discretizations of Radiative-Conductive-Diffusive Model

As stated in Section 2, with the assumptions that the thermal energy of the atmosphere is dominated by the gas phase and haze particles are in a thermal equilibrium with the gas, we can simplify the radiative-conductive-diffusive equation for the gas phase. Based on Equation (3) in the main text, we move the $ρ_g c_{p,g}$ to the right-hand side and remove the subscript g:

$$ \frac{\partial T}{\partial t} = R_{\text{rad}}^{*} + R_{\text{cond}}^{*} + R_{\text{adb}}^{*} + R_{\text{eddy}}^{*}. \quad (A1) $$

Here we define the radiative term $R_{\text{rad}}^*$ by combining both the gas and haze contribution:

$$ R_{\text{rad}}^{*} = \frac{1}{\rho c_p} (R_{\text{rad,g}}(T) + R_{\text{rad,p}}(T)), \quad (A2) $$

where T, ρ, and c_p are the temperature, gas density, and gas specific heat, respectively.

The conduction term of temperature $R_{\text{cond}}^*$ in a spherical coordinate is:

$$ R_{\text{cond}}^{*} = \frac{1}{\rho c_p r^2} \frac{\partial}{\partial r} \left( r^2 K_c \frac{\partial T}{\partial r} \right), \quad (A3) $$

where r is the distance from the center of Pluto and K_c is the conductivity.

The adiabatic cooling term due to the escape of the atmosphere $R_{\text{adb}}^*$ is:

$$ R_{\text{adb}}^{*} = -\frac{1}{\rho c_p} \sum_{i} \frac{m_i F_i}{4\pi r^2} \left( c_{p,i} \frac{\partial T}{\partial r} + \frac{GM}{r^2} \right), \quad (A4) $$

where G is the gravitational constant, M is the mass of Pluto, and m_i and c_{p,i} are the molecular mass and specific heat of the escaping species i, respectively. F_i is the modified Jeans escape flux $F_i = \Gamma(\lambda_i)[4\pi r^2 N_i \sqrt{\frac{k_b T_{\text{exo}}}{2\pi m_i}} (1 + \lambda_i)e^{-\lambda_i}]_{r=r_{\text{exo}}}$, where λ_i is the Jeans parameter, N_i is the molecular number density, k_b is the Boltzmann constant, and T_{\text{exo}} is the temperature at the exobase $r_{\text{exo}}$. Γ(λ_i) is the enhancement factor for the particle thermal escape (Zhu et al. 2014).

Finally, the eddy-diffusion term of potential temperature $R_{\text{eddy}}^*$ is:

$$ R_{\text{eddy}}^{*} = \frac{1}{\rho c_p r^2} p^{\kappa} \frac{\partial}{\partial r} \left( r^2 \rho c_p K_{zz} p^{-\kappa} \left( \frac{\partial T}{\partial r} + \frac{g}{c_p} \right) \right), \quad (A5) $$

where K_{zz} is the eddy heat diffusivity, p is the pressure, and g is the gravity.

For the vertical discretization, we use an unequal step size following Strobel et al. (1996). The atmosphere from the ground to ~1800 km is divided into 176 uneven layers, and the layer thickness increases with height by a factor of 1.03 between adjacent layers. We define $\Delta r_j = (r_{j+1} - r_{j-1})/2$, $\Delta r_j^+ = r_{j+1} - r_j$, and $\Delta r_j^- = r_j - r_{j-1}$.

At each vertical level j, we discretize Equation (A1):

$$ \frac{T_j^{n+1} - T_j^n}{\Delta t} = \left[ R_{\text{rad}}^* \right]_j^n + \left[ R_{\text{cond}}^* + R_{\text{adb}}^* + R_{\text{eddy}}^* \right]_j^{n+\frac{1}{2}}, \quad (A6) $$

where $\Delta t = t^{n+1} - t^n$ is the time step. n + 1/2 represents the average between the time step n (explicit) and n + 1 (implicit) in the Crank-Nicolson method. The radiative heating and cooling term $R_{\text{rad}}^*$ including multiple-scattering radiative transfer calculation needs to be treated explicitly. Other individual heat transfer terms on the right-hand side are treated as follows.

### A.1. Heat Conduction

We define K = r²K_c(T), K_{j+1/2} = $\frac{K_{j+1} + K_j}{2}$, and M_j = (ρ_j c_p r_j² Δr_j)⁻¹. We further assume $K_{j+1/2}^n \approx K_{j+1/2}^{n+1}$ and $M_j^n \approx M_j^{n+1}$. The heat conduction term $R_{\text{cond}}^*$ can be discretized as:

$$ \begin{aligned} \left[ R_{\text{cond}}^* \right]_j^{n+\frac{1}{2}} &= \frac{1}{2} \left( \left[ \frac{1}{\rho c_p r^2} \frac{\partial}{\partial r} \left( K \frac{\partial T}{\partial r} \right) \right]_j^{n+1} \right. \\ & \qquad \left. + \left[ \frac{1}{\rho c_p r^2} \frac{\partial}{\partial r} \left( K \frac{\partial T}{\partial r} \right) \right]_j^n \right) \\ &= \frac{M_j^n K_j^n}{2 \Delta r_j^+} (T_{j+1}^{n+1} - T_j^{n+1}) \\ & \qquad - \frac{M_j^n K_{j-1}^{n-1}}{2 \Delta r_j^-} (T_j^{n+1} - T_{j-1}^{n+1}) \\ & \qquad + \frac{M_j^n K_{j+1}^{n-1}}{2 \Delta r_j^+} (T_{j+1}^n - T_j^n) \\ & \qquad - \frac{M_j^n K_{j-1}^n}{2 \Delta r_j^-} (T_j^n - T_{j-1}^n). \end{aligned} \quad (A7) $$
## A.2. Adiabatic Cooling

We define $D_j = -\frac{1}{\rho_j c_p 4\pi r_j^2} (\sum_i m_i F_i c_{p,i})$ and $E_j = -\frac{1}{\rho_j c_p 4\pi r_j^2} (\sum_i m_i F_i)g$, and the adiabatic cooling term $R_{adb}^*$ can be discretized as:

$$ 
\begin{aligned}
\left[ R_{adb}^* \right]_j^{n+\frac{1}{2}} &= \left[ \left( -\frac{1}{\rho c_p} \frac{1}{4\pi r^2} \left( \sum_i m_i F_i c_{p,i} \right) \right) \frac{\partial T}{\partial r} \right. \\
&\qquad \left. + \left( -\frac{1}{\rho c_p} \frac{1}{4\pi r^2} \left( \sum_i m_i F_i \right) \right) g \right]_j^{n+\frac{1}{2}} \\
&= \frac{D_j^n}{4\Delta r_j} (T_{j+1}^{n+1} - T_j^{n+1}) \\
&\quad + \frac{D_j^n}{4\Delta r_j} (T_j^{n+1} - T_{j-1}^{n+1}) \\
&\quad + \frac{D_j^n}{4\Delta r_j} (T_{j+1}^n - T_j^n) \\
&\quad + \frac{D_j^n}{4\Delta r_j} (T_j^n - T_{j-1}^n) + E_j^n.
\end{aligned} 
\quad (A8) $$

## A.3. Eddy Heat Diffusion

We define $H = r^2 \rho c_p K_{zz} p^{-\kappa}$, $H_{j+\frac{1}{2}} = \frac{H_{j+1} + H_j}{2}$, $N_j = K_{zz}/(H_j \Delta r_j)$, and $L_{j+\frac{1}{2}} = g H_{j+\frac{1}{2}} / c_p$. We further assume $H_{j+\frac{1}{2}}^n \approx H_{j+\frac{1}{2}}^{n+1}$, and the eddy term $R_{eddy}^*$ can be discretized as:

$$ 
\begin{aligned}
\left[ R_{eddy}^* \right]_j^{n+\frac{1}{2}} &= \left[ \frac{K_{zz}}{H} \frac{\partial}{\partial r} \left( H \frac{\partial T}{\partial r} \right) \right]_j^{n+\frac{1}{2}} \\
&\quad + \left[ \frac{K_{zz}}{H} \frac{\partial}{\partial r} \left( H \frac{g}{c_p} \right) \right]_j^{n+\frac{1}{2}} \\
&= \frac{N_j^n H^n}{2 \Delta r_j^{n+1}} (T_{j+1}^{n+1} - T_j^{n+1}) \\
&\quad - \frac{N_j^n H^n}{2 \Delta r_j^{n+1}} (T_j^{n+1} - T_{j-1}^{n+1}) \\
&\quad + \frac{N_j^n H^n}{2 \Delta r_j^{n+1}} (T_{j+1}^n - T_j^n) \\
&\quad - \frac{N_j^n H^n}{2 \Delta r_j^{n+1}} (T_j^n - T_{j-1}^n) \\
&\quad + N_j^n (L_{j+\frac{1}{2}}^n - L_{j-\frac{1}{2}}^n).
\end{aligned} 
\quad (A9) $$

## A.4. Formulate the Tridiagonal Linear System

Now we can incorporate the above heat transfer terms (Equations (A7), (A8), and (A9)) into Equation (A6) and

formulate a linear system. We define the following coefficients at time step $n$:

$$ 
\begin{aligned}
A_j &= \left( \frac{M_j^n K_{j+\frac{1}{2}}^n}{2 \Delta r_j^+} + \frac{D_j^n}{4 \Delta r_j} + \frac{N_j^n H_{j+\frac{1}{2}}^n}{2 \Delta r_j^-} \right) \Delta t, \\
B_j &= \left( \frac{M_j^n K_{j-\frac{1}{2}}^n}{2 \Delta r_j^-} - \frac{D_j^n}{4 \Delta r_j} + \frac{N_j^n H_{j-\frac{1}{2}}^n}{2 \Delta r_j^-} \right) \Delta t, \\
C_j &= \left( \frac{N_j^n K_{j+\frac{1}{2}}^n}{2 \Delta r_j^+} + \frac{D_j^n}{4 \Delta r_j} + \frac{N_j^n H_{j+\frac{1}{2}}^n}{2 \Delta r_j^+} (T_{j+1}^n - T_j^n) \right. \\
&\qquad \left. - \left( \frac{N_j^n K_{j-\frac{1}{2}}^n}{2 \Delta r_j^-} - \frac{D_j^n}{4 \Delta r_j} + \frac{N_j^n H_{j-\frac{1}{2}}^n}{2 \Delta r_j^-} (T_j^n - T_{j-1}^n) \right) \Delta t \right. \\
&\qquad \left. + (R_{rad}^* + E_j^n + N_j^n (L_{j+\frac{1}{2}}^n - L_{j-\frac{1}{2}}^n)) \Delta t \right).
\end{aligned} 
\quad (A10) $$

Then the linear system is:

$$ 
\begin{aligned}
&-B_j T_{j-1}^{n+1} + (1 + A_j + B_j) T_j^{n+1} \\
&-A_j T_{j+1}^{n+1} = C_j + T_j^n.
\end{aligned} 
\quad (A11) $$

To solve this system, we let $a_j = -B_j$, $b_j = (1 + A_j + B_j)$, and $c_j = -A_j$, $g_j = C_j + T_j^n$. The system can be written using a tridiagonal matrix:

$$ 
\left[ \begin{array}{cccccccccc}
b_1 & c_1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\
a_2 & b_2 & c_2 & 0 & \cdots & 0 & 0 & 0 & 0 \\
0 & a_3 & b_3 & c_3 & \cdots & 0 & 0 & 0 & 0 \\
\cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
\cdots & \cdots & \cdots & \cdots & \cdots & a_{N-2} & b_{N-2} & c_{N-2} & 0 \\
0 & 0 & 0 & \cdots & 0 & a_{N-1} & b_{N-1} & c_{N-1} & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 & a_N & b_N & 0
\end{array} \right] \\
\times \left[ \begin{array}{c} T_1 \\ T_2 \\ T_3 \\ \cdots \\ T_{N-2} \\ T_{N-1} \\ T_N \end{array} \right] = \left[ \begin{array}{c} g_1 \\ g_2 \\ g_3 \\ \cdots \\ g_{N-2} \\ g_{N-1} \\ g_N \end{array} \right]^n.
$$

The boundary conditions are fixed bottom temperature and no heat flux at the top:

$$ 
\left\{ 
\begin{array}{l}
T_{r=r_0} = T_s \Rightarrow T_1 = T_s \Rightarrow b_1 = 1, c_1 = 0, g_1 = T_s, \quad (A12) \\
\frac{\partial T_\infty}{\partial r} = 0 \Rightarrow T_N = T_{N-1} \Rightarrow a_N = 1, b_N = -1, g_N = 0, \quad (A13)
\end{array}
\right. $$
such that

$$ \begin{bmatrix} 1 & 0 & 0 & 0 & \dots & \dots & 0 & 0 & 0 \\ a_2 & b_2 & c_2 & 0 & \dots & \dots & 0 & 0 & 0 \\ 0 & a_3 & b_3 & c_3 & \dots & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & \dots & a_{N-2} & b_{N-2} & c_{N-2} & 0 \\ 0 & 0 & 0 & \dots & \dots & 0 & a_{N-1} & b_{N-1} & c_{N-1} \\ 0 & 0 & 0 & \dots & \dots & 0 & 0 & 1 & -1 \end{bmatrix}_W $$

$$ \times \begin{bmatrix} T_1 \\ T_2 \\ T_3 \\ \vdots \\ T_{N-2} \\ T_{N-1} \\ T_N \end{bmatrix}^{n+1} = \begin{bmatrix} T_s \\ g_2 \\ g_3 \\ \vdots \\ g_{N-2} \\ g_{N-1} \\ 0 \end{bmatrix}^n . $$

Finally, the tridiagonal matrix (W) can be inverted numerically (e.g., using the Gaussian elimination method) to calculate the temperature profile of Pluto.

$$ \begin{bmatrix} T_1 \\ T_2 \\ T_3 \\ \vdots \\ T_{N-2} \\ T_{N-1} \\ T_N \end{bmatrix}^{n+1} = W^{-1} \begin{bmatrix} T_s \\ g_2 \\ g_3 \\ \vdots \\ g_{N-2} \\ g_{N-1} \\ 0 \end{bmatrix}^n . $$

References

Anderson, C. M., Samuelson, R. E., & Nna-Mvondo, D. 2018, SSRv, 214, 125

Bertrand, T., & Forget, F. 2016, Natur, 540, 86

Bertrand, T., & Forget, F. 2017, Icar, 287, 72

Bertrand, T., Forget, F., Schmitt, B., White, O. L., & Grundy, W. M. 2020a, NatCo, 11, 5056

Bertrand, T., Forget, F., White, O., et al. 2020b, JGRE, 125, e06120

Cheng, A. F., Summers, M. E., Gladstone, G. R., et al. 2017, Icar, 290, 112

Colaïtis, A., Spiga, A., Hourdin, F., et al. 2013, JGRE, 118, 1468

Etling, D., & Brown, R. A. 1993, BoLMe, 65, 215

Forget, F., Bertrand, T., Vangvichith, M., et al. 2017, Icar, 287, 54

Forget, F., Hourdin, F., Fournier, R., et al. 1999, JGR, 104, 24155

Gao, P., Fan, S., Wong, M. L., et al. 2017, Icar, 287, 116

Gladstone, G. R., Stern, S. A., Enrico, K., et al. 2016, Sci, 351, aad8866

Gladstone, G. R., & Young, L. A. 2019, AREPS, 47, 119

Hillier, J., Helfenstein, P., Verbiscer, A., & Veverka, J. 1991, JGRS, 96, 19203

Hinson, D. P., Linscott, I. R., Young, L. A., et al. 2017, Icar, 290, 96

Hudson, R. L., Gerakines, P. A., & Moore, M. H. 2014, Icar, 243, 148

Jacobs, A. D., Summers, M. E., Cheng, A. F., et al. 2021, Icar, 356, 113825

Jovanović, L., Gautier, T., Broch, L., et al. 2021, Icar, 362, 114398

Jovanović, L., Gautier, T., Vuitton, V., et al. 2020, Icar, 346, 113774

Khanna, R. K., Ospina, M. J., & Zhao, G. 1988, Icar, 73, 527

Khare, B. N., Sagan, C., Arakawa, E. T., et al. 1984, Icar, 60, 127

Krasnopolsky, V. A. 2020, Icar, 335, 113374

Krasnopolsky, V. A., & Cruikshank, D. P. 1999, JGR, 104, 21979

Lavvas, P., Lellouch, E., Strobel, D. F., et al. 2021, NatAs, 5, 289

Lavvas, P., Yelle, R. V., Koskinen, T., et al. 2013, PNAS, 110, 2729

Lellouch, E., Gurwell, M., Butler, B., et al. 2017, Icar, 286, 289

Luspay-Kuti, A., Mandt, K., Jessup, K.-L., et al. 2017, MNRAS, 472, 104

Mandt, K., Luspay-Kuti, A., Hamel, M., et al. 2017, MNRAS, 472, 118

Martonchik, J. V., & Orton, G. S. 1994, ApOpt, 33, 8306

Moore, M. H., Ferrante, R. F., Moore, W. J., & Hudson, R. 2010, ApJS, 191, 96

Ramirez, S. I., Coll, P., da Silva, A., et al. 2002, Icar, 156, 515

Stern, S. A., Bagenal, F., Enrico, K., et al. 2015, Sci, 350, aad1815

Stern, S. A., Kammer, J. A., Barth, E. L., et al. 2017, AJ, 154, 43

Strobel, D. F., & Zhu, X. 2017, Icar, 291, 55

Strobel, D. F., Zhu, X., Summers, M. E., & Stevens, M. H. 1996, Icar, 120, 266

Wong, M. L., Fan, S., Gao, P., et al. 2017, Icar, 287, 110

Yelle, R. V., Lunine, J. I., & Hunten, D. M. 1991, Icar, 89, 347

Young, L. A., Kammer, J. A., Steffl, A. J., et al. 2018, Icar, 300, 174

Zalucha, A. M., Zhu, X., Gulbis, A. A. S., Strobel, D. F., & Elliot, J. L. 2011, Icar, 214, 685

Zhang, X., & Showman, A. P. 2018, ApJ, 866, 1

Zhang, X., Strobel, D. F., & Imanaka, H. 2017, Natur, 551, 352

Zhu, X., Strobel, D. F., & Erwin, J. T. 2014, Icar, 228, 301