# Disk-resolved Photometric Properties of Pluto and the Coloring Materials across its Surface

Silvia Protopapa, Cathy B. Olkin, Will M. Grundy, Jian-Yang Li, Anne Verbiscer, Dale P. Cruikshank, Thomas Gautier, Eric Quirico, Jason C. Cook,

Dennis Reuter, et al.

► To cite this version:

Silvia Protopapa, Cathy B. Olkin, Will M. Grundy, Jian-Yang Li, Anne Verbiscer, et al. Disk-resolved Photometric Properties of Pluto and the Coloring Materials across its Surface. The Astronomical Journal, 2020, 159 (2), pp.74. 10.3847/1538-3881/ab5e82. insu-02468367

HAL Id: insu-02468367

https://insu.hal.science/insu-02468367v1

Submitted on 3 Aug 2025

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

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

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

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

# Disk-resolved Photometric Properties of Pluto and the Coloring Materials across its Surface

Silvia Protopapa<fnref n="1" />, Cathy B. Olkin<fnref n="1" />, Will M. Grundy<fnref n="2" />, Jian-Yang Li<fnref n="3" />, Anne Verbiscer<fnref n="4" />, Dale P. Cruikshank<fnref n="5" />, Thomas Gautier<fnref n="6" />, Eric Quirico<fnref n="7" />, Jason C. Cook<fnref n="8" />, Dennis Reuter<fnref n="9" />, Carly J. A. Howett<fnref n="1" />, Alan Stern<fnref n="1" />, Ross A. Beyer<fnref n="5" />,<fnref n="10" />, Simon Porter<fnref n="1" />, Leslie A. Young<fnref n="1" />, Hal A. Weaver<fnref n="1" />, Kim Ennico<fnref n="5" />, Cristina M. Dalle Ore<fnref n="5" />,<fnref n="10" />, Francesca Scipioni<fnref n="10" />, and Kelsi Singer<fnref n="1" />

*Received 2019 August 27; revised 2019 November 29; accepted 2019 December 1; published 2020 January 31*

## Abstract

A multiwavelength regionally dependent photometric analysis of Pluto's anti-Charon-facing hemisphere using images collected by *New Horizons*' Multispectral Visible Imaging Camera (MVIC) reveals large variations in the absolute value and spectral slope of the single-scattering albedo. Four regions of interest are analyzed: the dark equatorial belt, Pluto's north pole, nitrogen-rich regions, and the mid-latitude terrains. Regions dominated by volatile ices such as Lowell Regio and Sputnik Planitia present single-scattering albedos of ~0.98 at 492 nm, almost neutral across MVIC's visible wavelength range (400–910 nm), indicating limited contributions from tholin materials. Pluto's dark equatorial regions, informally named Cthulhu and Krun Maculae, have single-scattering albedos of ~0.16 at 492 nm and are the reddest regions. Applying the Hapke radiative transfer model to combined MVIC and Linear Etalon Imaging Spectral Array (LEISA) spectra (400–2500 nm) of Cthulhu Macula and Lowell Regio successfully reproduces the spectral properties of these two regions of dramatically disparate coloration, composition, and morphology. Since this model uses only a single coloring agent, very similar to the Titan-like tholin of Khare et al., to account for all of Pluto's colors, this result supports the Grundy et al. conclusion that Pluto's coloration is the result of photochemical products mostly produced in the atmosphere. Although cosmic rays and extreme ultraviolet photons reach Pluto's surface where they can drive chemical processing, observations of diverse surface colors do not require different chemical products produced in different environments. We report a correction scaling factor in the LEISA radiometric calibration of 0.74 ± 0.05.

*Unified Astronomy Thesaurus concepts: Trans-Neptunian objects (1705); Pluto (1267); Radiative transfer simulations (1967)*

*Supporting material: data behind figures*

## 1. Introduction

The chemistry of Pluto's atmosphere and surface has become a key factor in understanding the origin and evolution of this icy dwarf planet, and by extension that of a vast number of similar sized and smaller bodies in the Kuiper Belt, beyond the terrestrial and giant planets. A major role in the composition of the surface and atmosphere of Pluto is played by macro-molecular organic material analogous to laboratory tholins (from the Greek word meaning "sepia ink"). These are the residues produced through the action of several energy sources on small organic molecules (e.g., Khare et al. 1984; Cruikshank 2005).

An ample variety of tholins can be synthesized in the laboratory from native materials found in Pluto's atmosphere (Gladstone et al. 2016) and surface (Grundy et al. 2016a; Protopapa et al. 2017). Macromolecular carbonaceous solids relevant to solar system bodies have been produced in the laboratory for many years using CH₄, N₂, and other simple

molecular starting recipes, in the gas phase, exposed to electrical discharge in cold plasma conditions (e.g., Khare et al. 1984; Szopa et al. 2006; Quirico et al. 2008; Imanaka et al. 2012), and an ultraviolet (UV) emission line lamp (e.g., Tran et al. 2003; Vuitton et al. 2009). These tholins are representative of gas phase chemistry and aerosols formation. In contrast, UV photons or charged-particle irradiation experiments on icy materials (e.g., Materese et al. 2015) are more representative of the solid-state chemistry of the surface exposed to UV photons and cosmic rays. Overall, laboratory tholins absorb in the UV-visible and near-infrared (NIR) spectral range and present different degrees of coloration in the visible wavelength range (see Figure 6 by Mahjoub et al. 2012 for an example of tholins color variation). The spectral behavior of a tholin depends on the recipe used in the laboratory to generate the material (e.g., composition and phase of the initial mixture, energy source, pressure).

In enhanced color images, Pluto's vivid red, brown, and yellow colors paint complex patterns that appear to be

<footnote n="1">¹Southwest Research Institute Boulder, CO 80302, USA; sprotopapa@boulder.swri.edu</footnote>

<footnote n="2">²Lowell Observatory Flagstaff, AZ 86001, USA</footnote>

<footnote n="3">³Planetary Science Institute Tucson, AZ 85719, USA</footnote>

<footnote n="4">⁴University of Virginia Charlottesville, VA 22904, USA</footnote>

<footnote n="5">⁵NASA Ames Research Center Moffett Field, CA 94035, USA</footnote>

<footnote n="6">⁶LATMOS-IPSL, CNRS, Sorbonne Université, UVSQ Guyancourt F-78280, France</footnote>

<footnote n="7">⁷Université Grenoble Alpes, CNRS, Institut de Planetologie et Astrophysique de Grenoble (IPAG) UMR 5274 Grenoble, F-38041, France</footnote>

<footnote n="8">⁸Pinhead Institute Telluride, CO, USA</footnote>

<footnote n="9">⁹NAŞA Goddard Space Flight Center Greenbelt, MD, USA</footnote>

<footnote n="10">¹⁰SETI Institute Mountain View, CA 94043, USA</footnote>

<footnote n="11">¹¹Johns Hopkins University Applied Physics Laboratory Laurel, MD 20723, USA</footnote>
correlated with Pluto's highly varied underlying geological structures (Stern et al. 2015; Olkin et al. 2017). The color contrast is less obvious in natural-color images. A consensus has emerged around the concept that tholins are present on the surface of Pluto and serve as coloring agents (e.g., Stern et al. 2018). However, the number of distinct types of tholins on the surface of Pluto, and the processes responsible for their formation and distribution remain the subject of investigation.

Olkin et al. (2017) showed the existence of two distinct color trends in visible wavelength color–color space, strongly favoring the existence of at least two distinct colorants. Grundy et al. (2018) argued that dark red photochemical haze settles out of the atmosphere globally, so it ought to be possible to account for all of Pluto's colors with a single pigment. However, they were unable to match Pluto's diverse colors with models based on a single pigment, concluding that if photochemical haze particles were indeed responsible for all of Pluto's coloration, the particles must be further modified in certain environments subsequent to settling out of the atmosphere.

We analyze the spectral signature(s) of Pluto's colored terrain(s) to characterize and intercompare the potentially different coloring agents and ultimately understand their history and chemical evolution. We use Pluto data from the *New Horizons* Ralph instrument (Reuter et al. 2008), which consists of a visible multicolor/panchromatic mapper, the Multispectral Visible Imaging Camera (MVIC), and a mapping IR composition spectrometer, the Linear Etalon Imaging Spectral Array (LEISA). Our investigation provides a systematic analysis of Pluto's non-icy components across the surface. This analysis helps understand the origin and nature of such materials in the local frame as well as with respect to the trans-Neptunian object population, which presents a variety of colorations reminiscent of those found on Pluto.

## 2. Pluto's Disk-resolved Photometric Properties in the Visible Wavelength Range

We report here a detailed study of disk-resolved photometric properties of Pluto using *New Horizons* Ralph/MVIC images in the visible wavelength range of 400–910 nm acquired during the Pluto flyby (the initial results and details of the flyby can be found in Stern et al. 2015). Regions across Pluto's closest approach hemisphere are investigated in this paper. The derivation of a regionally based photometric model permits us to (1) decouple the intrinsic surface albedo variability from effects related to the observing geometry, and therefore investigate quantitatively the true heterogeneity of Pluto's surface; (2) combine visible (MVIC) and NIR observations (LEISA) accounting for the different viewing geometries at which these data were acquired; and (3) model visible and NIR measurements of Pluto to derive quantitative information of its surface composition.

### 2.1. MVIC Pluto Data

MVIC is equipped with seven CCD arrays. Four of these enable color mapping in the visible wavelength range (Table 1). For details about the instrument, its performance, and data calibration, the reader is referred to Reuter et al. (2008), Howett et al. (2017), and Olkin et al. (2017).

<table><caption>**Table 1**
MVIC Filter Passbands</caption><thead><tr><td>Filter</td><td>Wavelength Range (nm)</td><td>Pivot wavelength (nm)</td></tr></thead><tbody><tr><td>Blue</td><td>400–550</td><td>492</td></tr><tr><td>Red</td><td>540–700</td><td>624</td></tr><tr><td>NIR</td><td>780–975</td><td>861</td></tr><tr><td>CH4</td><td>860–910</td><td>883</td></tr></tbody></table>

Note. The pivot wavelengths of the filters are from Howett et al. (2017).

The physical positions of the four MVIC color CCDs are offset from one another, so the target (due to scanning and spacecraft and target motion) comes into their view at slightly different times. Therefore, each color channel needs to be registered separately to the relevant Pluto base map. This registration was performed by Schenk et al. (2018). With this registration, the Integrated Software for Imagers and Spectroscopy (ISIS; US Geological Survey 2019) was used to orthographically project all MVIC color images to a common perspective view, which is that of a sphere at the target's size and location relative to the spacecraft as of the mid-scan header time of the blue filter image. The reprojection was done to a target grid with a spatial scale higher than the native MVIC pixel scale. The point of subsampling is to minimize degradation of spatial information due to the nearest neighbor resampling. The Pluto latitude and longitude as well as the viewing geometry (defined by the incident angle $i$, emergent angle $e$, and phase angle $g$) of each pixel of the MVIC color scan under consideration is computed. It is worth pointing out that we do not account for Pluto's local topography, given that the projection is done onto a sphere. This implies that small misregistrations between the color channels associated with steep topography can occur. The incidence and emergent angles do not account for local topography either.

While observations of Pluto were acquired with MVIC from 2015 April through July, we focus here on the five color images taken within the last two days on approach to Pluto (Table 2). Each line in the table describes a single MVIC color scan that covered Pluto. The request ID in Table 2 is an identifier of the observation used in planning the observing sequence while the mission elapsed time (MET) is a counter on the spacecraft that increments each second since launch and provides a unique timestamp for the data. The MET at the start of each MVIC scan is reported in Table 2.

The derivation of a regionally based photometric model requires the selection of the same surface regions across data acquired at different viewing geometries. This requirement calls for the five MVIC scans listed in Table 2 and displayed in Figure 1. Pluto's observations recorded prior to those selected, cover a different subspacecraft longitude due to Pluto's rotation (period of approximately 6.4 days) or were acquired at a spatial resolution too poor for our purposes.

MVIC images used for this analysis are calibrated to units of radiance factor (RADF). This is also commonly referred to as 1/F and it is defined as the ratio of the bidirectional reflectance of a surface $r$ to that of a perfectly diffuse surface of the same size and distance to the Sun and observer, but illuminated at a normal direction (Hapke 1993). Several radiometric keywords can be used to convert *New Horizons* MVIC data from number
<table><caption>**Table 2**
MVIC Observations of Pluto</caption><thead><tr><th>Request ID</th><th>MET</th><th>Mid-time of Observation UTC</th><th>Image Scale (km pixel<sup>-1</sup>)</th><th>Sub-s/c Long (Deg E)</th><th>Sub-s/c Lat (Deg N)</th></tr></thead><tbody><tr><td>PEMV_01_PCNH_Multi_Long_1d1</td><td>299104958</td><td>2015 Jul 13 T14:50:51.941</td><td>20.78</td><td>196.53</td><td>42.70</td></tr><tr><td>PEMV_01_PC_Multi_Long_1d2</td><td>299127628</td><td>2015 Jul 13 T21:08:40.941</td><td>14.53</td><td>182.07</td><td>42.50</td></tr><tr><td>PEMV_01_PC_Color_TimeRes</td><td>299147983</td><td>2015 Jul 14 T02:47:54.441</td><td>8.92</td><td>169.48</td><td>42.06</td></tr><tr><td>PEMV_01_PC：<br>Color_1</td><td>299162518</td><td>2015 Jul 14 T06:50:11.942</td><td>4.91</td><td>161.41</td><td>41.12</td></tr><tr><td>PEMV_01_P：<br>Color2</td><td>299178098</td><td>2015 Jul 14 T11:10:52.442</td><td>0.66</td><td>168.03</td><td>25.98</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Five Pluto images from *New Horizons*. The spatial resolution of the images improves with the image number because each successive image was taken at a closer range to Pluto than the previous one. These enhanced color maps were produced with MVIC's radiance factor (RADF) images in the blue, red, and NIR filter displayed in the blue, green, and red color channels, respectively. Regions where the incidence angle from the Sun or the emergent angle to the spacecraft exceeds 80° from the zenith were excluded in the analysis. The viewing geometry of each scan and spatial resolution are listed in Table 2.</figcaption></figure>

of counts (DN) to RADF. The radiometric keyword RSOLAR is adopted for our data set. For more information about the procedure applied to convert DN to RADF and the difference between the several keywords, the reader is referred to Howett et al. (2017).

We resize each MVIC scan to the dimensions of the first Pluto observation at the lowest spatial resolution (MET 299104958). This way all scans have comparable spatial

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Long Range Reconnaissance Imager (LORRI) panchromatic image of Pluto's surface with latitude and longitude markings and formal and informal feature names.</figcaption></figure>

resolution and the same region of interest (ROI) can be easily selected across the measurements acquired with different viewing geometries. Each color image in a given scan is convolved with a Gaussian kernel with a FWHM equal to that of each Ralph MVIC channel—the average of the in-track and cross-track FWHM of the MVIC color channels is computed, see Reuter et al. (2008).

The next step of our analysis is the determination of the error on the RADF. The derived error will be used in the best-fit minimization of the RADF data to determine Pluto's photometric properties. For each pixel, the standard deviation of its neighboring pixels located in a 3-by-3 pixel box is computed and assumed to be the 1σ RADF error. The 3σ error is used in the best-fit minimization processes.

## 2.2. Regions of Interest

For our study of disk-resolved photometric properties, we have identified four ROIs across Pluto's encounter hemisphere (see Figure 2): (a) the terrains of Cthulhu<fnref n="12" /> and Krun Maculae, (b) the yellow hue on Pluto's north pole, (c) the volatile-rich terrains including, but not limited to, Sputnik Planitia, and (d)

<footnote n="12">¹² Some names in use on Pluto and Charon are now formalized and others are still informal.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Panel A: enhanced color image of Pluto (PEMV_01_PCNH_Multi_Long_1d1_01) with MVIC's blue, red, and NIR filter images displayed in the blue, green, and red color channels, respectively. Panel B: the same as panel A, with selected ROIs shown with four different colors (red, yellow, magenta, and green; see the text for details). Panel C: distribution of NIR/red and red/blue color ratios. The Sun symbol indicates neutral colors; redder colors extend up and to the right. The selected pixels shown in panel B and representative of the four ROIs are displayed in the color-color diagram with the same color convention.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Panel (a) shows the photometric MVIC data in the blue channel for ROI a as a function of viewing geometry. The top, middle, and bottom panels display the RADF as a function of phase angle, incidence angle, and emergent angle, respectively. Colored symbols are used to identify data from each MVIC scan. Panels (b) through (d) are the same as panel (a) but for ROIs b through d, respectively.</figcaption></figure>

the less red terrains of Baret Montes and Viking Terra to the north of Cthulhu and south of 35°N. These regions are shown in Figure 3 (panel B) overlaid with different colors on the

PEMV_01_PCNH_Multi_Long_1d1 color image. Red, yellow, magenta, and green points represent ROIs a, b, c, and d, respectively. The ROI selection resulted from the analysis of
the Pluto color variability by Olkin et al. (2017). Regions that cluster in the color–color diagram and therefore present similar color properties have been selected. As evident in panel C of Figure 3, the four ROIs span almost entirely the full color–color diagram. The areas of Cthulhu and Krun Maculae (ROI a, red points) and the volatile-rich terrains (Protopapa et al. 2017; Schmitt et al. 2017b), like Sputnik Planitia (ROI c, magenta points), correspond to the end members on the top right side and bottom left side of the steep-sloped mixing line of the color–color diagram, respectively, while the regions of Baret Montes and Viking Terra (ROI d, green points) fall among the intermediate pixels along this mixing line. The yellow hue on Pluto’s north pole (ROI b, yellow points) lays across the less-steep color mixing line of the color–color diagram. Analysis of these four ROIs will provide a sense of the range of Pluto’s photometric properties in the visible wavelength range. The data we fit models to are shown in Figure 4, which displays the reduced photometric data in the blue filter as a function of viewing geometry for each ROI. Data points with i > 80° or e > 80° are discarded to avoid pixels too close to the limb or terminator.

## 2.3. Photometric Modeling

Photometric models describe the dependence of the RADF on scattering geometry, which is defined by the incident angle $i$, emergent angle $e$, and phase angle $g$. The RADF is equal to $\pi r$, where $r$ is the surface bidirectional reflectance. We adopt the Hapke multiple-scattering radiative transfer model (Hapke 2002, 2012) to derive the photometric properties of Pluto. This is the same model used by Protopapa et al. (2017) to investigate the spatial distribution of the abundance and textural properties of the volatile and nonvolatile materials across Pluto’s surface. For self-consistency of this paper we briefly summarize the modeling below. The RADF is given by

$$ \text{RADF}(i, e, g) = \frac{w}{4} \frac{\mu_{0e}}{\mu_{0e} + \mu_e} \{[1 + B(B_0, h, g)]p(\xi, g) + H(\mu_{0e}, w)H(\mu_e, w) - 1\}S(i, e, g, \theta), \quad (1) $$

where $w$ is the single-scattering albedo, defined as the ratio of the scattering efficiency to the extinction efficiency. Therefore, $w$ ranges from 0 to 1, where $w = 0$ implies that the particles absorb all the radiation. $S$ is the shadowing function and accounts for effects due to macroscopic roughness $\theta$, and $\mu_{0e}$ and $\mu_e$ are the cosine of the incident and emergent angle, respectively, corrected for roughness.

The backscatter function $B$ accounts for the shadow-hiding opposition effect (the nonlinear increase in the reflectance phase curve with decreasing phase angle) and is given by

$$ B(g) = \frac{B_0}{1 + 1/h \tan(g/2)}, \quad (2) $$

where $B_0$ and $h$ represent the amplitude and width of the opposition effect, respectively. We do not account for the coherent backscattering opposition effect (CBOE) given the limited phase angle coverage of our data (Figure 4). CBOE generally appears at phase angles <2°, while our data have a minimum phase angle of ~14° (Hapke 2002).

$H(\mu, w)$ is the Chandrasekhar $H$-function and the term $H(\mu_{0e}, w)H(\mu_e, w) - 1$ in Equation (1) characterizes multiple

scattering assuming isotropic scatterers. We adopt the $H$-function approximation suggested by Hapke (2002).

Finally, $p(\xi, g)$ is the single-particle phase function and describes the angular pattern into which the radiation is scattered by the surface. A single lobe Henyey–Greenstein function (Henyey & Greenstein 1941) is used for $p(\xi, g)$ and is given by

$$ p(\xi, g) = \frac{1 - \xi^2}{(1 + 2\xi \cos g + \xi^2)^{3/2}}, \quad (3) $$

where $\xi$ is the cosine asymmetry factor, which can range from -1 to 1. The scattering is isotropic ($p(g) = 1$) when $\xi = 0$. If $\xi > 0$, $p(g)$ increases monotonically between 0 and $\pi$, and decreases monotonically if $\xi < 0$.

We use the Levenberg–Marquardt $\chi^2$ minimization algorithm to find the best-fit solutions (Markwardt 2009). The model quality is quantified by means of the reduced $\chi^2$, defined as

$$ \chi_{\nu}^{2} = \frac{1}{\nu} \sum_{i=0}^{N} \left( \frac{\text{RADF}_{i} - \text{RADF}_{i, \text{model}}(B_{0}, h, w, \xi, \theta)}{\text{ERR}_{i}} \right)^{2}, \quad (4) $$

where RADF$_i$, ERR$_i$, and RADF$_{i, model}$ represent the measured RADF, the error on the RADF, and the modeled RADF, respectively, in correspondence of the individual pixel on the surface identified by the subscript $i$. The modeling depends from $B_0$, $h$, $w$, $\xi$, and $\theta$, while the geometry ($i, e, g$) of each pixel is known. The sum is over all $N$ data points. The degree of freedom, $\nu$, equals the number of observations $N$ minus the number of fitted parameters $m$.

## 2.4. Results

We conduct a multiwavelength photometric analysis of Pluto, meaning that each MVIC color filter is treated separately. This is because $w$, $\xi$, $B_0$, and $h$ are wavelength dependent. On the contrary, $\theta$, being related to surface morphology, has no spectral dependence.

Given the limited phase angle coverage, we consider the amplitude of the opposition effect $B_0$ as well as the compaction parameters $h$ constant throughout our wavelength range and across the entire surface of Pluto. We adopt values of 0.307 and 0.206 for $B_0$ and $h$, respectively, obtained from a disk-integrated analysis of the Hubble Space Telescope data at $\lambda = 0.55\mu$m (Verbiser et al. 2019). We solve instead for single-scattering albedo $w$, cosine asymmetry factor $\xi$, and mean roughness slope $\theta$.

Figure 5 shows the joint error distribution in the blue filter of ($w, \xi$), ($w, \theta$), and ($\xi, \theta$) for the dark terrains of Cthulhu Macula and Krun Macula. These maps were obtained computing the reduced $\chi^2$ over a grid of parameter values, where the grid spans a parameter range much wider than the eventual errors. For each point in the grid, which defines the values of the two parameters under consideration, we solve for the third parameter, with $B_0$ and $h$ still fixed. The maps in Figure 5 clearly show that the fits are not sensitive to the mean roughness slope $\theta$, in part due to the limited phase angle coverage. Therefore, we set $\theta = 20°$ (Verbiser et al. 2019), and let $w$ and $\xi$ free. Our approach is very similar to that used by Hillier et al. (1994) to assess, using Voyager 2 data, the photometric diversity of terrains on Triton.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 2D reduced $\chi^2$ distributions resulting from the parameter exploration of the fit to the data points in ROI a (Figure 3, panel B, red points) in the blue filter. The solid curves shown in white, red, and yellow with $\Delta\chi^2 = 2.30, 4.61, 9.21$ represent the 68%, 90%, and 99% confidence levels, respectively. We conclude that the phase angle range is too limited to solve for the mean roughness slope $\theta$.</figcaption></figure>

The regions of Cthulhu and Krun Maculae are considered among the darkest terrains of Pluto's encounter hemisphere. This is based on the RADF being very low between 0.05 and 0.2 (see Figure 4, panel a). However, looking at the RADF values alone is not sufficient to support this statement. As an example, the volatile-rich terrains on Pluto (ROI c) display RADF as low as 0.2 at incident angles of 80° (see Figure 4(c), middle panel). This observation is concrete evidence that the intrinsic surface albedo variability needs to be decoupled from effects related to the observing geometry.

To ensure that the solution obtained from the best-fit minimization algorithm does not correspond to a local minimum rather than the global minimum and investigate the possible intercorrelation between the Hapke parameters, we explore the full parameter space $(w, \xi)$ with $w$ ranging between 0 and 1 with a step of 0.01 and $\xi$ between -1 and 1 with a step of 0.01. We determine the reduced $\chi^2$ distributions and the degree of confidence in the best-fit values. The joint error distribution of $(w, \xi)$ in all filters for ROIs a, b, and c is shown in Figure 6. The 68%, 90%, and 99% confidence region contours are displayed in white, red, and yellow, respectively. White lines/squares show the best-fit parameters. Their values and the correspondent $1\sigma$ errors are listed in Table 3 and shown in Figure 7. There is evidence for a strong positive correlation between $\xi$ and $w$, but both parameters are well constrained, such that we can safely rely on the wavelength trend and absolute value of both $w$ and $\xi$.

The goodness of fit for ROI c is plotted in Figure 8. By comparing the reduced $\chi^2$ among different ROIs in Table 3, it is clear that ROI c represents the worst case scenario in terms of best-fit modeling. However, we attribute the large reduced $\chi^2$ to the wide range of geologic terrains encompassed by this ROI

including the flat ice sheet of Sputnik Planitia as well as the bladed terrains along the low latitude eastern terminator limb of Pluto's encounter hemisphere (Moore et al. 2018). In spite the relatively high reduced-$\chi^2$, the measured RADF and the RADF at the same observation geometry modeled with the derived Hapke parameters (RADF<sub>model</sub>) are linearly correlated, as expected. This result is nominal and it is the same for ROI a and b as well. We point out that considering a less extended region with respect to ROI c yields similar results, but with larger error bars.

ROIs a, b, and c are backscattering, with $\xi$ constant, within the error, throughout the wavelength range, and equal to an average of $-0.21 \pm 0.07$. Given the small number of points in ROI d, we set $\xi = -0.21$ fixed and solve for $w$ (see Figure 7, green points). Notice that our data were acquired at phase angles ranging between ~15° and ~40° (see Figure 4), where the backscattering lobe dominates. Therefore, it is not unexpected to find $\xi$ to be negative from this analysis. Nevertheless, it is worth noticing that our value for $\xi$ is relatively close to the global value of $-0.36$ at $0.55~\mu\text{m}$ found by Verbischer et al. (2019) for Pluto and $-0.33$ at $0.56~\mu\text{m}$ by Hillier (1994) for Triton.

## 2.5. Discussion

The region including Cthulhu Macula and Krun Macula is a very distinctive surface unit (ROI a, Figure 7), being significantly darker and redder than the other three. Less red and less dark terrains are distributed at mid-latitudes (e.g., Baret Montes and Viking Terra, ROI d). The polar cap terrains (ROI b) as well as the volatile-rich terrains including but not limited to Sputnik Planitia (ROI c) appear neutral in color since they display a single-scattering albedo almost flat with wavelengths. Also, they are extremely bright with single-particle scattering albedos of ~0.98 or higher. Notice that the yellow coloration of Pluto's north pole is due to the RADF in the blue filter being suppressed with respect to that in the red, NIR, and CH4. This trend is also visible in the spectral behavior of the single-scattering albedo when accounting for the errors.

Given that tholin materials absorb strongly in the visible wavelength range (e.g., Brassé et al. 2015; low albedo and red slope), regions with $w$ as high as 0.98 or higher must be highly dominated by ice materials and depleted of coloring agents. Therefore, we attribute variations in albedo and spectral slope between different terrains across the surface of Pluto to changes in abundance and grain size of these organic materials and not necessarily to different coloring agents. This idea is supported by the tholin abundance map presented by Protopapa et al. (2017) and shown in Figure 9: a high concentration of dark compounds is found in Cthulhu Macula while Lowell Regio is highly depleted in tholins. The composition maps by Protopapa et al. (2017) are the result of a pixel-by-pixel Hapke radiative transfer model applied to *New Horizons* LEISA scans. Because important compositional information exists in the visible wavelength range where low albedo organic compounds present the most diagnostic spectral signatures, the true contribution of the coloring agents cannot be assessed if the visible spectral domain is disregarded. Therefore, our idea that a single color pigment accounts for all of Pluto's colors in different concentrations and particle size requires a radiative transfer model of the visible and NIR measurements of Pluto.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Joint error distribution of $(w, \xi)$ in all filters for ROI a, b, and c is shown in the top, middle, and bottom panel, respectively. The solid curves shown in white, red, and yellow with $\Delta\chi^2 = 2.30, 4.61, 9.21$ represent the 68%, 90%, and 99% confidence levels, respectively. White lines/squares show the best-fit parameters.</figcaption></figure>

The next section presents the results of this analysis for Cthulhu Macula and Lowell Regio.

**3. Tholin Characterization across Pluto: Two Provinces**

**3.1. MVIC and LEISA Radiometric Calibration, Charon as a Calibrator**

Our analysis requires the use of LEISA and MVIC co-registered data, in order to span the full wavelength range between 0.4 and 2.5 µm, where both ices and tholins present diagnostic spectral features and behaviors. As a preliminary

step in the analysis we investigate the radiometric calibration of MVIC and LEISA, using Charon as a test case.

Pluto's largest satellite Charon appears fairly homogeneous in NIR scans recorded by the Ralph/LEISA instrument on board of the *New Horizons* spacecraft (Grundy et al. 2016a; Dalle Ore et al. 2018). LEISA spectroscopic measurements in the NIR wavelength range display strong water-ice and ammonia features, and an overall continuum that becomes blue at wavelengths longward of about ~1.4 µm. This spectrum is consistent with ground-based observations reported by several authors (e.g., Buie & Grundy 2000; Dumas et al. 2001; Cook et al. 2007;
<table><caption>**Table 3**
Best-fit Parameters for Pluto's ROIs (See the Text for Details) with $B_0 = 0.307$, $h = 0.206$, $\theta = 20^\circ$ Fixed in the Model Fitting</caption><thead><tr><th>ROI</th><th>Pivot Wavelength (nm)</th><th>w</th><th>ξ</th><th>χ<sup>2</sup><sub>ν</sub></th></tr></thead><tbody><tr><td rowspan="4">ROI a</td><td>492</td><td>0.16<sup>+0.07</sup><sub>-0.04</sub></td><td>-0.36<sup>+0.16</sup><sub>-0.10</sub></td><td>2.16</td></tr><tr><td>624</td><td>0.32<sup>+0.14</sup><sub>-0.08</sub></td><td>-0.29<sup>+0.18</sup><sub>-0.10</sub></td><td>1.93</td></tr><tr><td>861</td><td>0.61<sup>+0.11</sup><sub>-0.19</sub></td><td>-0.21<sup>+0.13</sup><sub>-0.09</sub></td><td>3.69</td></tr><tr><td>883</td><td>0.67<sup>+0.11</sup><sub>-0.10</sub></td><td>-0.17<sup>+0.15</sup><sub>-0.10</sub></td><td>2.78</td></tr><tr><td rowspan="4">ROI b</td><td>492</td><td>0.98<sup>+0.01</sup><sub>-0.09</sub></td><td>-0.07<sup>+0.19</sup><sub>-0.24</sub></td><td>0.95</td></tr><tr><td>624</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.22<sup>+0.02</sup><sub>-0.06</sub></td><td>1.59</td></tr><tr><td>861</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.21<sup>+0.01</sup><sub>-0.06</sub></td><td>2.30</td></tr><tr><td>883</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.17<sup>+0.02</sup><sub>-0.07</sub></td><td>1.92</td></tr><tr><td rowspan="4">ROI c</td><td>492</td><td>0.97<sup>+0.02</sup><sub>-0.03</sub></td><td>-0.22<sup>+0.16</sup><sub>-0.09</sub></td><td>6.27</td></tr><tr><td>624</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.24<sup>+0.03</sup><sub>-0.06</sub></td><td>8.02</td></tr><tr><td>861</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.20<sup>+0.02</sup><sub>-0.07</sub></td><td>13.06</td></tr><tr><td>883</td><td>0.99<sup>+0.00</sup><sub>-0.01</sub></td><td>-0.14<sup>+0.02</sup><sub>-0.09</sub></td><td>7.69</td></tr><tr><td rowspan="4">ROI d</td><td>492</td><td>0.62<sup>+0.08</sup><sub>-0.08</sub></td><td>-0.21<sup>a</sup></td><td>0.28</td></tr><tr><td>624</td><td>0.75<sup>+0.05</sup><sub>-0.06</sub></td><td>-0.21<sup>a</sup></td><td>0.45</td></tr><tr><td>861</td><td>0.85<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.21<sup>a</sup></td><td>1.45</td></tr><tr><td>883</td><td>0.85<sup>+0.01</sup><sub>-0.02</sub></td><td>-0.21<sup>a</sup></td><td>1.51</td></tr></tbody></table>

**Note.**

Verbiser et al. 2007). Charon's surface is spectrally neutral in MVIC color ratios, with the exception of the northern polar cap characterized by red coloration possibly associated to the presence of tholins (Grundy et al. 2016a, 2016b). Ground-based visible spectroscopy of Charon, separated from Pluto, reveals a completely featureless and almost perfectly flat spectrum across the wavelength range of 0.5–1 µm (Fink & Disanti 1988). While there is no published spectrum of Charon across the full wavelength range of 0.4–2.5 µm, the evidence reported above suggest that, with the exception of the red pole, Charon displays an overall flat continuum, turning blue beyond ~1.4 µm. This fairly simple spectral behavior is ideal to verify the radiometric calibration of the visible and IR components of the Ralph instrument.

We consider MVIC and LEISA scans of Charon acquired consecutively. As a preliminary step in the analysis, for each MVIC and LEISA pair, we reproject each MVIC color frame to the perspective view of the target as seen from the spacecraft at the mid-scan header time of the LEISA scan. Notice that each set of MVIC and LEISA measurements share almost the same viewing geometry and no corrections for incident, emergent, and phase angle needs to be applied when combining the two data sets (see Figure 10). The characteristics of each MVIC and LEISA pair used for this analysis are reported in Table 4. The viewing geometry (i, e, g) as well as latitude and longitude reported in Table 4 are the average across the full visible disk, after the reprojection is applied. A representative MVIC and LEISA combined I/F spectrum of Charon's surface, outside of the red polar cap, is shown in Figure 11 (top panel). Contrary to the expectation, an offset between the continuum in the visible and IR wavelength ranges is observed. This offset is quantified by performing the ratio between the average continuum level computed in the range of 1.37–1.39 µm and in the MVIC NIR and narrow band CH4 visible channels. This factor is fairly constant across all the MVIC+LEISA pairs (Figure 11, bottom panel), which supports the idea of a systematic error in the

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Single-scattering albedo *w* and cosine asymmetry factor *ξ* as a function of wavelengths for four ROIs across Pluto's encounter hemisphere.</figcaption></figure>

MVIC or LEISA radiometric calibration. We adopt the average value of 1.35 ± 0.09 as a scaling to adjust the LEISA measurements to match the continuum level in the visible (e.g., LEISA spectra need to be divided by 1.35 ± 0.09 or multiplied by 0.74 ± 0.05.). We apply the correction factor to the LEISA data, rather than to the MVIC data, given the accurate in-flight radiometric calibration of MVIC observations. This calibration was based on two semi-independent techniques, one of which used Charon as a calibrator. For details see Howett et al. (2017).

## 3.2. Pluto: Combined MVIC and LEISA Measurements

Given Pluto's disk-resolved photometric properties, it is now possible to combine MVIC (visible range 0.40–0.91 µm) and LEISA observations (NIR range 1.25–2.50 µm) accounting for the different viewing geometries at which these data were acquired. This combination is done by computing, for each ROI, the visible photometrically corrected radiance factor RADF', which is the measured RADF at $i_{VIS}$, $e_{VIS}$, $g_{VIS}$, reported to the LEISA geometry ($i_{IR}$, $e_{IR}$, $g_{IR}$) by means of the

<footnote>a This parameter has been set as constant.</footnote>
<table><caption>**Table 4**
MVIC and LEISA Observations of Charon</caption><thead><tr><th>Instrument</th><th>Request ID</th><th>MET</th><th>Mid-time of Observation (UTC)</th><th>g (deg)</th><th>e (deg)</th><th>i (deg)</th><th>Lat (deg)</th><th>Lon (deg)</th></tr></thead><tbody><tr><td>MVIC</td><td>PC_MULTI_MAP_B_15 MC</td><td>299025872</td><td>2015 Jul 12 T16:53:04.940</td><td>15.2</td><td>44.5</td><td>46.1</td><td>38.8</td><td>113.8</td></tr><tr><td>LEISA</td><td>PC_MULTI_MAP_B_15 LE</td><td>299026199</td><td>2015 Jul 12 T17:01:49.374</td><td>15.2</td><td>45.9</td><td>47.4</td><td>38.4</td><td>113.2</td></tr><tr><td>MVIC</td><td>PC_MULTI_MAP_B_17 MC</td><td>299064592</td><td>2015 Jul 13 T03:38:06.440</td><td>15.5</td><td>44.3</td><td>45.7</td><td>38.9</td><td>130.1</td></tr><tr><td>LEISA</td><td>PC_MULTI_MAP_B_17 LE</td><td>299064869</td><td>2015 Jul 13 T03:46:31.874</td><td>15.5</td><td>45.6</td><td>47.2</td><td>38.9</td><td>127.8</td></tr><tr><td>MVIC</td><td>PC_MULTI_MAP_B_18 MC</td><td>299079022</td><td>2015 Jul 13 T07:38:36.941</td><td>15.6</td><td>44.5</td><td>44.1</td><td>39.3</td><td>157.2</td></tr><tr><td>LEISA</td><td>PC_MULTI_MAP_B_18 LE</td><td>299079314</td><td>2015 Jul 13 T07:46:25.375</td><td>15.6</td><td>45.8</td><td>45.2</td><td>39.0</td><td>158.5</td></tr><tr><td>MVIC</td><td>PCNH_Multi_Long_1d01 MC</td><td>299104952</td><td>2015 Jul 13 T14:50:51.941</td><td>16.1</td><td>44.6</td><td>46.3</td><td>38.8</td><td>158.4</td></tr><tr><td>LEISA</td><td>PCNH_Multi_Long_1d1 LE</td><td>299105209</td><td>2015 Jul 13 T14:59:32.875</td><td>16.1</td><td>44.9</td><td>46.7</td><td>38.8</td><td>158.1</td></tr><tr><td>MVIC</td><td>PC_Color-TimeRes</td><td>299147977</td><td>2015 Jul 14 T02:47:54.441</td><td>18.3</td><td>44.7</td><td>47.0</td><td>38.5</td><td>194.2</td></tr><tr><td>LEISA</td><td>C_Leisa</td><td>299146219</td><td>2015 Jul 14 T02:21:49.875</td><td>18.1</td><td>44.9</td><td>47.1</td><td>38.5</td><td>193.3</td></tr><tr><td>MVIC</td><td>PC_Color_1</td><td>299162512</td><td>2015 Jul 14 T06:50:11.942</td><td>21.1</td><td>45.0</td><td>48.6</td><td>36.9</td><td>209.0</td></tr><tr><td>LEISA</td><td>C_LEISA_LORRI_1 LE</td><td>299171308</td><td>2015 Jul 14 T09:20:28.376</td><td>27.0</td><td>45.1</td><td>49.9</td><td>36.9</td><td>208.7</td></tr><tr><td>MVIC</td><td>C_Color_2</td><td>299176432</td><td>2015 Jul 14 T10:42:28.942</td><td>38.8</td><td>45.3</td><td>53.7</td><td>35.0</td><td>204.8</td></tr><tr><td>LEISA</td><td>C_LEISA_Hires LE</td><td>299175509</td><td>2015 Jul 14 T10:32:05.876</td><td>35.8</td><td>45.2</td><td>53.6</td><td>35.0</td><td>204.6</td></tr></tbody></table>

Note. g: phase angle; e: emergent angle; i: incident angle; Lat: latitude; and Lon: longitude. The values reported for i, e, g, Lon, and Lat are the average across the full visible disk.

<figure><img src="image_10.png" /><figcaption>**Figure 11.** Top panel: Charon spectrum over the wavelength range of 0.4–2.5 µm. The spectrum has been obtained by averaging the reflectance spectra in the blue box of Figure 10. The 1σ standard deviation is shown. Bottom panel: NIR to visible ratio computed for Charon spectra extracted out of each MVIC+LEISA pair considered in this analysis. The horizontal solid line represents the mean of the seven values, while the dashed lines show the 1σ standard deviation.</figcaption></figure>

Cook et al. (2019) as “Cthulhu Regio H₂O-rich,” but corrected to account for the radiometric calibration.

### 3.3. Modeling

We use the scattering radiative transfer model of Hapke (2012), described in Section 2.3. We consider an areal or an intimate mixture, or a combination of the two. In the case of areal mixture, the bidirectional reflectance spectra of the individual surface components ($r_j$) are summed with weights equal to the fractional area of each terrain ($F_j$). In the case of an intimate mixture the averaging process is at the level of the individual particle. To model an intimate mixture of different compositions, we need to modify the parameters of the particles that enter in the RADF equation (Equation (1), Section 2.3). As an example, the volume single-scattering albedo becomes

$$ \bar{w} = \left[ \sum_j \frac{V_j}{D_j} w_j \right] \left[ \sum_j \frac{V_j}{D_j} \right]^{-1}, \qquad (6) $$

where the subscript j refers to any property of the particles (e.g., composition), and $V_j$ and $D_j$ are the fractional volume and particle diameter of the jth type of particle, respectively. Therefore, the parameters that enter in the reflectance equation are averages of the properties of the various types of particles in the mixture weighted by the cross-sectional area. We adopt the equivalent slab model presented by Hapke (2012) to compute $w_j$, given the complex refractive index and the path length of the material. This approach is valid in the geometric-optics regime when the particles are much larger than the wavelengths ($\frac{\pi D}{\lambda} > 1$). For small particles, we follow the approach adopted by Cook et al. (2019), which in turn relies on the formulation described by Clark et al. (2010). The Hapke parameters listed in Table 3 are adopted to model the spectra of Cthulhu (ROI a) and Lowell Regio (ROI b), respectively. A constant value of
<figure><img src="image_11.png" /><figcaption>**Figure 12.** Left panel: top panel shows combined MVIC and LEISA mean RADF spectrum of eastern Cthulhu (red filled circles) compared with the best-fit model obtained considering two types of tholins, designated as sample A (the same as Titan tholin by Khare et al. (1984) but with *k* modified in the NIR, see the text for details) and sample B (the same as sample A but modified in the visible wavelength range; data and model have been offset along the y-axis by 0.2 for clarity). The models using sample A and B are shown as the dashed-dotted blue line and solid blue line, respectively. Open green diamonds and squares show the convolution of the best-fit models with the MVIC transmission filters. The bottom panel is the same as top panel but for the yellow hue on Pluto's north pole. Only the model obtained with sample B tholin is shown. Right panel: the single-scattering albedo of eastern Cthulhu (top panel, red filled circles) and Lowell Regio (bottom panel, yellow filled diamonds) is compared with the single-scattering predicted by the models (open green diamonds and squares). The MVIC and LEISA mean RADF spectra of eastern Cthulhu and Lowell Regio behind this figure are available as machine-readable tables.</figcaption></figure>

(The data used to create this figure are available.)

−0.17 in the NIR is adopted for the cosine asymmetry parameter $\xi$ in both cases. We assume that $h$, $B_0$, and $\theta$ are constant across the full wavelength range. The free parameters in our model are the effective diameter ($D_j$) and the contribution of each surface terrain to the mixture ($F_j, V_j$).

We model Cthulhu Macula considering an areal mixture of two surface units. The first (Unit 1) is an intimate mixture of tholin and aggregates of water ice with methane ice inclusions. The second is composed of tholin. The optical constants of the aggregate ($n(\lambda), k(\lambda)$) are estimated from effective medium theory (Bruggeman mixing formula; Bohren & Huffman 1983), where the effective medium is a mixture of water ice and methane. The use of methane ice imbedded into water ice enables the attenuation of the absorption bands of methane, which is otherwise spectroscopically dominant in the NIR, and

obtain a better fit to the data around 1.8 and 2.3–2.4 µm. The model we propose here does not require CH$_3$OH (Cook et al. 2019). While methane provides a satisfactory fit and has been used also by Cook et al. (2019) to model the same region of Pluto, other hydrocarbons (Clark et al. 2009) could be alternative surface materials. Due to the lack of optical constants and because it is beyond the scope of this paper, we did not explore alternative hydrocarbons. Therefore, this model may not be a unique solution and depends on the surface components considered. It is worth noting that our main goal is to reproduce the overall spectral behavior and absolute values of MVIC+LEISA RADF measurements of Pluto. While weak absorptions are important, we are not focusing in this work on interpreting their nature. The spectral characteristic of the tholin material employed in our modeling are detailed below in this section.
<figure><img src="image_12.png" /><figcaption>**Figure 13.** Imaginary part of the refractive index, *k*, of each surface material in the wavelength range of 0.3–2.4 µm. The MVIC transmission filters are shown for reference. For details see the text. The optical constants of tholins sample A and B behind this figure are available as machine-readable tables.

(The data used to create this figure are available.)</figcaption></figure>

A challenge in modeling spectroscopic measurements throughout the visible and NIR resides in the lack of optical constants of analog materials in the full wavelength range. We describe below the optical constants used for each surface material (see Figure 13).

The most complete set of optical constants of crystalline water ice spanning the full range from the UV to the far-IR has been made by Warren & Brandt (2008). However, most of the data used for the calculation of the optical constants have been recorded at fairly high temperatures (between 80 and 263 K) and are thus not fully relevant for trans-Neptunian object studies (de Bergh et al. 2008). Absorption coefficient spectra of hexagonal water ice at temperatures between 20 and 270 K in the range between 1.0 and 2.7 µm have been published by Grundy & Schmitt (1998). We compute the imaginary part of the refractive index *k* from the absorption coefficient spectrum *α* at 60 K reported by Grundy & Schmitt (1998) between 1.0 and 2.7 µm (using *α* = 4π*k*/λ) and combined it with that by Warren & Brandt (2008) below 1 µm. We use the real part of the refractive index *n* by Warren & Brandt (2008) in the full range of wavelengths.

Tabulated values of (*n*, *k*) for crystalline CH₄–I at 39 K (Grundy et al. 2002) in the range of 0.7–5.0 µm are available in

the Solid Spectroscopy Hosting Architecture of Databases and Expertise (SSHADE) database (Schmitt et al. 2017a).<fnref n="13" /> We set *k* equal to 0 in the range of 0.4–0.7 µm. This way we are assuming that no methane absorptions occur below 0.7 µm and we avoid introducing any artificial slopes. We extrapolate *n* to shorter wavelengths assuming a constant value equal to that tabulated at 0.7 µm. The optical constants of CH₄ in the solid solution in β-N₂ at 36.5 K (Quirico & Schmitt 1997) are available in the range of 1.0–5.0 µm. Given that MVIC data are not sensitive to the methane dilution state and therefore to the wavelength shift of the CH₄ absorption bands, we adopt the *k* values of pure CH₄ for CH₄ diluted in N₂ below 1.0 µm. The real part *n* of the refractive index is constant throughout the wavelength range and equal to the visible value of N₂ (*n* = 1.23, Quirico et al. 1999).

Titan and Triton tholins, which are produced in the laboratory from gaseous methane–nitrogen mixtures exposed to electrical discharge or energetic photons, are among the non-ice materials most widely used in modeling the colors and spectra of trans-Neptunian objects and Centaurs. Spectrally, the major difference between Triton tholin and Titan tholin is in the NIR, with the former being redder. Laboratory reflectance measurements have been obtained for refractory tholins particularly relevant to Pluto (Pluto's ice tholin) through UV and low-energy electron bombardment of a ice mixture representative of Pluto's average surface composition (N₂:CH₄:CO = 100:1:1; Materese et al. 2015; Cruikshank et al. 2016). The spectrum shows a red slope between 0.5 and 1 µm, and turns blue in the range between 1 and 2.5 µm. A preliminary set of optical constants have been obtained assuming a particle size for the tholin material, and used in the literature to model Pluto and Charon data in the NIR (Dalle Ore et al. 2018; Cook et al. 2019). We acknowledge that the size range can play a crucial role in the derivation of the imaginary part of the refractive index *k*. In spite of Pluto's ice tholin being a potential good analog to model Pluto's data (this has not yet been demonstrated given that Pluto's ice tholin was applied to model NIR data without accounting for the visible component of Pluto's spectrum), we prefer here to rely on published values of (*n*, *k*) and we thus disregard Pluto's ice tholin. We adopt instead Titan atmospheric tholins in our modeling given the neutral slope of the Pluto LEISA spectrum of Cthulhu Macula.

Many studies have been performed to improve the experimental database of the optical constants of Titan atmospheric tholins (Brassé et al. 2015). Most of the studies characterize the refractive indices in the visible wavelength range (Ramirez et al. 2002; Imanaka et al. 2004; Mahjoub et al. 2012; Sciamma-O'Brien et al. 2012) or in the mid-IR range (Imanaka et al. 2012). However, few publications report the determination of the complex refractive indices of Titan tholins in the wavelength range of 0.4–2.5 µm except for Khare et al. (1984) and Tran et al. (2003). As observed by Brassé et al. (2015), the *k* values reported by Khare et al. (1984) and Tran et al. (2003) present the same global trend of *k* decreasing with increasing wavelengths from 0.20 to 2.5 µm. However, Tran et al. (2003), who consider four different kinds of Titan haze analogs, report *k* values in the NIR (0.7–2.5 µm) that are two orders of magnitude higher with respect to those by Khare et al. (1984). The Khare et al. (1984) measurements are also higher than other laboratory tholins at other wavelengths (Brassé et al.

<footnote n="13">¹³ https://www.sshade.eu/</footnote>
<table><caption>**Table 5**
Models</caption><thead><tr><th rowspan="3">Region</th><th colspan="8">Unit 1</th><th colspan="2">Unit 2</th></tr><tr><th colspan="2">H<sub>2</sub>O + CH<sub>4</sub></th><th colspan="2">Sample B</th><th colspan="2">CH<sub>4</sub>:N<sub>2</sub></th><th colspan="2">N<sub>2</sub>:CH<sub>4</sub></th><th rowspan="2">F<sub>N<sub>2</sub>:CH<sub>4</sub></sub> (%)</th><th rowspan="2">D(µm)</th><th rowspan="2">F<sub>unit2</sub>(%)</th></tr><tr><th>V(%)</th><th>D(µm)</th><th>V<sup>CH<sub>4</sub></sup><sub>H<sub>2</sub>O+CH<sub>4</sub></sub> (%)</th><th>V(%)</th><th>D(µm)</th><th>V(%)</th><th>D(µm)</th><th>V(%)</th><th>D(µm)</th></tr></thead><tbody><tr><td>Cthulhu Macula</td><td>0.2</td><td>0.1</td><td>89.7</td><td>99.8</td><td>4.6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.01</td><td>17.7</td></tr><tr><td>Lowell Regio</td><td>...</td><td>...</td><td>...</td><td>0.002</td><td>0.8</td><td>41.1</td><td>733.5</td><td>58.8</td><td>7280.3</td><td>0.15</td><td>...</td><td>...</td></tr></tbody></table>

2015). For a detailed discussion about the main causes behind such a discrepancy, the reader is referred to Brassé et al. (2015). Here we notice that the set of optical constants reported by Tran et al. (2003) present *k* values too high in the NIR with respect to that in the visible and are not adequate to reproduce Pluto measurements. Therefore, we use the (*n*, *k*) values for Titan tholin reported by Khare et al. (1984).

The *k* values of Titan tholin by Khare et al. (1984) present a sharp increase beyond 1.8 µm, which implies the presence of an absorption band centered around ~2.2 µm. However, such an absorption feature is not displayed by our Pluto data extracted in correspondence of Cthulhu Macula. We point out that this feature in the laboratory Titan tholin by Khare et al. (1984) is very broad and extends from 1.8 to 2.4 µm. The sampling of the Khare et al. (1984) optical constants is very poor and the sharp, and the abrupt change of *k* at 1.8 µm seems very anomalous and inconsistent with the data. We recognize that a 2.2 µm feature, possibly due to a C ≡ N stretching mode as suggested by Cruikshank et al. (1991), is visible in the tholin reflectance measurements presented by Materese et al. (2015) and some tholin samples by Imanaka et al. (2012; depending on the tholin pressure formation conditions). However, it is absent from other laboratory tholin spectra (Tran et al. 2003). We point out that Cook et al. (2018, 2019) and Dalle Ore et al. (2018) model the LEISA data of Pluto, Charon, and Pluto's small satellites adopting the set of optical constants derived from the Materese et al. (2015) ice tholin after removing from the *k* values any feature through a best-fit line.

The mismatch between the laboratory data by Khare et al. (1984) and the New Horizons observations prompt us to compute a synthetic set of (*n*, *k*) for a tholin material that enables us to fit the data. We modify the *k* values by Khare et al. (1984) beyond 1.8 µm assuming them to be constant with wavelengths up to 2.4 µm and equal to the mean value computed between 1.35 and 1.45 µm. No changes are applied to the real part of the refractive index *n*. We designate this set of (*n*, *k*) as sample A. We stress that we are aware that tholins present on the surface of Pluto might display features in the NIR, which are actually absent from our sample A tholin, but here we focus on the photometric level of the tholin material. Inferring tholin's absorption bands in the NIR using modeling techniques of Pluto's data is very challenging given that the solution would strongly depend on the surface compounds assumed in the modeling. The model obtained using sample A as tholin material is shown in the top left panel of Figure 12 (dashed-dotted blue line). Open green diamonds show the convolution of the best-fit model with the MVIC transmission filters, which are displayed in Figure 13. The model fits the data satisfactorily across the full wavelength range but in the MVIC RED filter centered at 0.62 µm. The discrepancy between data and the model is also visible when comparing the single-scattering albedo of Cthulhu, computed when solving for the

disk-resolved photometric properties of Pluto (Section 2) and the modeled single-scattering albedo. However, the model is still within $1\sigma$ from the data (Figure 12, right top panel).

The effective optical constants (*k*) of tholin sample A are modified at 0.69 µm to better match the data, leading to tholin sample B. Notice that we preserve the original sampling of the data by Khare et al. (1984); this is why we modify *k* at 0.69 µm. The *k* values of tholins sample A and B compared to the original Khare et al. (1984) measurements are shown in Figure 13 (bottom panel). A lower *k* value for the tholin material in the red filter with respect to that tabulated by Khare et al. (1984) is needed to match the data (see Figure 12, top panel, solid blue line and open green squares). This ad hoc set of optical constants for tholin sample B enables us to fit the data not only in correspondence of Cthulhu Macula but also those of Lowell Regio. The comparison between the data of the yellow hue on Pluto's north pole and the best-fit model, obtained considering an intimate mixture of tholin (sample B) and N₂:CH₄ and CH₄:N₂ (see Protopapa et al. 2017 for details about N₂:CH₄ and CH₄:N₂), is shown in the bottom panel of Figure 12. The details of the best-fit modeling of Cthulhu Macula and Lowell Regio are listed in Table 5.

## 4. Conclusions and Discussion

We present a detailed modeling of the photometric properties of four ROIs across Pluto's surface at four visible wavelengths from 400 to 910 nm using the color images collected by the Ralph/MVIC instrument on board of New Horizons. With phase angle coverage limited to phase angles between ~14° and ~40°, we were unable to solve for the amplitude B₀ and width *h* of the opposition effect, as well as the roughness parameter *θ*, so we had to fix their values to those determined by Verbiser et al. (2019) of 0.307, 0.206, 20°, respectively. We solved for the single-scattering albedo *w* and the cosine asymmetry parameter *ξ*. The best-fit Hapke model parameters at four wavelengths together with their formal statistical uncertainties are listed in Table 3. Strong variations are observed in the absolute value (between ~0.16 and ~0.98 at 0.5 µm) and spectral slope of the single-scattering albedo. Cthulhu and Krun Maculae are significantly darker and redder than the rest of the surface. Regions dominated by volatile ices such as Sputnik Planitia and the yellow material across Pluto's north pole observed in enhanced color images present single-scattering albedos of ~0.98 or higher, and almost neutral across the visible wavelength range. This result indicates a very limited contribution of tholin materials on the optically active surfaces in these regions.

Through analysis of Ralph/MVIC and Ralph/LEISA data of Charon, we report that LEISA measurements need to be scaled down by a constant factor of 1.35 ± 0.09. The absolute level of RADF is a strong constraint when applying radiative transfer models to derive abundance and grain size composition
information. Therefore, previous quantitative analyses of LEISA data not accounting for this scaling factor should be revisited (e.g., Protopapa et al. 2017; Cook et al. 2019).

We analyze combined MVIC and LEISA mean RADF spectra of eastern Cthulhu and Lowell Regio (yellow material on Pluto's north pole). We model, for each region, both the RADF in the full wavelength range of 0.4-2.4 µm and the single-scattering albedo in the four MVIC filters, accounting for the photometric properties discussed above. Modeling not only the RADF but also *w* provides an additional constraint on the composition. While the RADF depends on the geometry, the photometric properties and the composition (abundance and grain size of each material) of the surface, the single scattering albedo depends exclusively on the composition of the surface. We use a common tholin material to reproduce the spectral properties of these two regions that have diverse colorations and compositions. With respect to several coloring agents for which optical constants are tabulated in the literature, the Titan tholin by Khare et al. (1984) provides a better match to the photometric level (absolute values of RADF) and spectral properties of Pluto's terrains. However, while the tholin material used to reproduce Pluto's spectral properties shares many common points with the *n*, *k* of Khare et al. (1984), two main differences appear: (1) the imaginary part of the refractive index that we find for Pluto's tholin is constant beyond 1.8 µm and does not display a peak at 2.2 µm and (2) a slight depression in the *k* value at 0.69 µm better fits the observations. Because a single pigment can be used to account for all of Pluto's colors, we concur with the idea suggested first by Grundy et al. (2018) that Pluto's coloration is the result of photochemical products mostly produced in the atmosphere, that condense as haze particles, and eventually settle to the surface. We consider less likely instead that tholins are produced through in situ radiolysis or photolysis of local ices, when the surface is not shielded from UV photons and solar wind particles. In fact, the striking diversity of Pluto's composition would respond differently to irradiation, as shown by a broad range of laboratory measurements (e.g., Brassé et al. 2015) and coloring agents with different spectral properties would be produced. If our hypothesis is correct and the color pigment originates as haze, then the various Titan tholins might be more appropriate analog materials, since they are produced from gas phase precursors. This work demonstrates the need for laboratory measurements reproducing the overall trend and absolute level of the real and imaginary part of the refractive index of our ad hoc tholin material to support our hypothesis.

In this paper we investigated large-scale units with obvious differences in tholin abundances (e.g., Cthulhu Macula and Lowell Regio). The same technique can be applied to smaller scale geologic regions such as the convection cells within Sputnik Planitia and the layered mountain blocks in the chaos region of the al-Idrisi mountains to investigate the age dependence and age variations of the atmospheric precipitation of colorants. As an example, it is apparent at high spatial resolution that there is a color difference between the center and the boundary of Sputnik Planum convection cells. Because convection theory (e.g., McKinnon et al. 2016) predicts a younger to older age gradient over several $10^5$ yr from cell center to margins, it is possible that future work could establish a correlation that would provide new support for the colorant precipitation hypothesis. Similarly, the layered colorant in the al-Idrisi mountain blocks provide evidence for varying colorant

deposition rates with epoch. These and similar new forms of related evidence will be investigated in future work.

This work was funded by NASA, through the New Frontiers Data Analysis program, via contract 80NSSC19K0821 to SwRI and by the *New Horizons* project. The authors gratefully acknowledge the *New Horizons* team and the support of NASA's *New Horizons* project for providing very high quality calibrated Ralph/MVIC and LEISA data, archived at the NASA's Planetary Data System Small Bodies Node (PDS-SBN). S.P. thanks the NASA grant 80NSSC19K0402 for partial funding that supported her work. We thank the anonymous reviewer for comments that improved the paper.

*Software:* USGS Integrated Software for Imagers and Spectrometers (ISIS), Solid Spectroscopy Hosting Architecture of Databases and Expertise (SSHADE, Schmitt et al. 2017a).

## ORCID iDs

Silvia Protopapa https://orcid.org/0000-0001-8541-8550

Jian-Yang Li https://orcid.org/0000-0003-3841-9977

Anne Verbiscer https://orcid.org/0000-0002-3323-9304

Simon Porter https://orcid.org/0000-0003-0333-6055

Hal A. Weaver https://orcid.org/0000-0003-0951-7762

Kim Ennico https://orcid.org/0000-0002-8847-8492

Kelsi Singer https://orcid.org/0000-0003-3045-8445

## References

Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (Wiley: New York)

Brassé, C., Muñoz, O., Coll, P., et al. 2015, P&SS, 109, 159

Buie, M. W., & Grundy, W. M. 2000, Icar, 148, 324

Clark, R. N., Curchin, J. M., Hoefen, T. M., et al. 2009, JGRE, 114, E03001

Clark, R. N., Pieters, C. M., Taylor, L. A., et al. 2010, LPI, 41, 2337

Cook, J. C., Dalle Ore, C. M., Protopapa, S., et al. 2019, Icar, 331, 148

Cook, J. C., Desch, S. J., Roush, T. L., et al. 2007, ApJ, 663, 1406

Cook, J. C., Ore, C. M. D., Protopapa, S., et al. 2018, Icar, 315, 30

Cruikshank, D. P. 2005, SSRv, 116, 421

Cruikshank, D. P., Allamandola, L. J., Hartmann, W. K., et al. 1991, Icar, 94, 345

Cruikshank, D. P., Stern, S. A., Grundy, W. M., et al. 2016, LPI, 47, 1696

Dalle Ore, C. M., Protopapa, S., Cook, J. C., et al. 2018, Icar, 300, 21

de Bergh, C., Schmitt, B., Moroz, L. V., et al. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 483

Dumas, C., Terrile, R. J., Brown, R. H., et al. 2001, AJ, 121, 1163

Fink, U., & Disanti, M. A. 1988, AJ, 95, 229

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

Grundy, W. M., Bertrand, T., Binzel, R. P., et al. 2018, Icar, 314, 232

Grundy, W. M., Binzel, R. P., Buratti, B. J., et al. 2016a, Sci, 351, aad9189

Grundy, W. M., Cruikshank, D. P., Gladstone, G. R., et al. 2016b, Natur, 539, 65

Grundy, W. M., & Schmitt, B. 1998, JGR, 103, 25809

Grundy, W. M., Schmitt, B., & Quirico, E. 2002, Icar, 155, 486

Hapke, B. 1993, Topics in Remote Sensing, 1st Edition (Cambridge: Cambridge Univ. Press)

Hapke, B. 2002, Icar, 157, 523

Hapke, B. 2012, Topics in Remote Sensing, 2nd Edition (Cambridge: Cambridge Univ. Press)

Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70

Hillier, J., Veverka, J., Helfenstein, P., et al. 1994, Icar, 109, 296

Howett, C. J. A., Parker, A. H., Olkin, C. B., et al. 2017, Icar, 287, 140

Imanaka, H., Cruikshank, D. P., Khare, B. N., et al. 2012, Icar, 218, 247

Imanaka, H., Khare, B. N., Elsila, J. E., et al. 2004, Icar, 168, 344

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

Lorenzi, V., Pinilla-Alonso, N., Licandro, J., et al. 2016, A&A, 585, A131

Mahjoub, A., Carrasco, N., Yahoo, P.-R., et al. 2012, Icar, 221, 670

Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler (San Francisco, CA: ASP), 251
Materese, C. K., Cruikshank, D. P., Sandford, S. A., et al. 2015, ApJ, 812, 150

McKinnon, W. B., Nimmo, F., Wong, T., et al. 2016, Natur, 534, 82

Moore, J. M., Howard, A. D., Umurhan, O. M., et al. 2018, Icar, 300, 129

Olkin, C. B., Spencer, J. R., Grundy, W. M., et al. 2017, AJ, 154, 258

Protopapa, S., Grundy, W. M., Reuter, D. C., et al. 2017, Icar, 287, 218

Quirico, E., Douté, S., Schmitt, B., et al. 1999, Icar, 139, 159

Quirico, E., Montagnac, G., Lees, V., et al. 2008, Icar, 198, 218

Quirico, E., & Schmitt, B. 1997, Icar, 127, 354

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

Reuter, D. C., Stern, S. A., Scherrer, J., et al. 2008, SSRv, 140, 129

Schenk, P. M., Beyer, R. A., McKinnon, W. B., et al. 2018, Icar, 314, 400

Schmitt, B., Bollard, P., Albert, D., et al. 2017a, SSHADE: “Solid Spectroscopy Hosting Architecture of Databases and Expertise,” OSUG Data Center, SSHADE (OSUG Data Center), Service/Database, doi:10.26302/SSHADE

Schmitt, B., Philippe, S., Grundy, W. M., et al. 2017b, Icar, 287, 229

Sciamma-O’Brien, E., Dahoo, P.-R., Hadamcik, E., et al. 2012, Icar, 218, 356

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

Stern, S. A., Grundy, W. M., McKinnon, W. B., et al. 2018, ARA&A, 56, 357

Szopa, C., Cernogora, G., Boufendi, L., et al. 2006, P&SS, 54, 394

Tran, B. N., Joseph, J. C., Ferris, J. P., et al. 2003, Icar, 165, 379

US Geological Survey 2019, Integrated Software for Imagers and Spectrometers (ISIS), https://www.usgs.gov/software/integrated-software-imagers-and-spectrometers-isis

Verbiser, A., Showalter, M., Helfenstein, P., et al. 2019, EPSC Abstracts, 13, EPSC-DPS2019-1261

Verbiser, A. J., Peterson, D. E., Skrutskie, M. F., et al. 2007, LPI, 38, 2318

Vuitton, V., Tran, B. N., Persans, P. D., et al. 2009, Icar, 203, 663

Warren, S. G., & Brandt, R. E. 2008, JGRD, 113, D14220