Testing tholins as analogues of the dark reddish material covering Pluto’s Cthulhu region

Abstract Pluto’s fly-by by the New Horizons spacecraft in July 2015 has revealed a dark reddish equatorial region, named Cthulhu, covered by a dark, non-icy material whose origin and composition have yet to be determined. It has been suggested that this material could form from the sedimentation of photochemical aerosols, originating from dissociation and ionisation processes in Pluto’s high atmosphere (similarly to aerosols forming Titan’s haze). This hypothesis is here further investigated by comparing New Horizons spectra collected both in the visible and the near-infrared to laboratory reflectance measurements of analogues of Pluto’s aerosols (Pluto tholins). These aerosols were synthesised in conditions mimicking Pluto’s atmosphere, and their optical and reflectance properties were determined, before being used in Hapke models. In particular, the single scattering albedo and phase function of Pluto tholins were retrieved through Hapke model inversion, performed from laboratory reflectance spectra collected under various geometries. From reconstructed reflectance spectra and direct comparison with New Horizons data, some of these tholins are shown to reproduce the photometric level (i.e. reflectance continuum) reasonably well in the near-infrared. Nevertheless, a misfit of the red visible slope still remains and tholins absorption bands present in the modelled spectra are absent in those collected by the New Horizons instruments. Several hypotheses are considered to explain the absence of these absorption features in LEISA data, namely high porosity effects or GCR irradiation. The formation of highly porous structures, which is currently our preferred scenario, could be promoted by either sublimation of ices initially mixed with the aerosols, or gentle deposition under Pluto’s weak gravity.


Introduction
Pluto is the largest of the Trans-Neptunian objects (TNOs). Before the flyby performed by the New Horizons' spacecraft in July, 2015, knowledge of Pluto's surface composition was limited to groundbased and Hubble Space Telescope (HST), disc-integrated observations. However, inferring Pluto's surface composition is of primary importance to gain insight into its bulk composition and to understand the interactions between the surface and Pluto's tenuous atmosphere. Before the New Horizons mission, CH 4 (Cruikshank et al., 1976), N 2 and CO ices (Owen et al., 1993) were identified from ground-based In addition to bright terrains covered by volatile ices, dark regions were already detected by ground-based and HST observations (e.g. Buie and Tholen, 1989;Buie et al., 2010;Grundy and Buie, 2002) and were later confirmed by the New Horizons' flyby (Grundy et al., 2016;Olkin et al., 2017;Protopapa et al., 2020). These regions are extremely low reflectance terrains showing a steep red slope in the visible. This advocates for the presence of a non-icy, dark material on Pluto's surface which has not yet been firmly identified. Most of these low reflectance terrains are dark red patchy regions located along Pluto's equatorial belt. The largest of these dark red patches, which extends from about 20 • east to 160 • east in longitude, is referred to as the Cthulhu region, and presents some of the reddest colours of Pluto's surface. Three main scenarios currently coexist to explain the presence of this dark material on the surface of Pluto.
First, this material could originate from the formation of the planet itself (Sekine et al., 2017). Pluto is thought to have formed from a giant impact that might have also created its satellite Charon (Canup, 2010;Desch, 2015;McKinnon et al., 2017). The extremely high energy delivered by such a collision between the proto-Pluto and the impactor is suspected to have triggered the melting of a significant part of the surface and icy-bedrock, creating a warm liquid water ocean. It has been shown that dark reddish, organic material can be synthesised in warm liquid water from simple organic compounds (e.g. formaldehyde) that are typically found in comets and could have been present on the proto-Pluto or brought in by the impactor (Sekine et al., 2017). Polymerisation of the simple organic compounds would lead to the formation of complex macromolecular compounds, responsible for the darkening and reddening of the water solution. Sample materials formed in warm water solution showed a steep red slope in the visible and no strong spectral feature in the near-infrared, which seems consistent with the New Horizons spectra collected over the dark red equatorial regions of Pluto (Sekine et al., 2017). Moreover, hydrodynamics simulation conducted in Sekine et al. (2017) shows that significant melted basins could have been created by a giant impact, the size of the largest one being comparable to that of the Cthulhu region.
Second, it has also been proposed that the red material could originate from the irradiation of ices on Pluto's surface by solar ultraviolet (UV) light and galactic cosmic rays (GCR), with a small contribution from charged solar wind particles (Cruikshank et al., 2015). As Pluto's tenuous atmosphere may have blocked most of the UV flux in the past, GCR was considered a more likely candidate, even if its energy flux at Pluto's distance is significantly lower than that of UV light. Laboratory experiments showed that both UV and charged particles irradiation of a mixture of the different ices present on Pluto (CH 4 , N 2 , CO, C 2 H 6 ) can trigger a rich chemistry and form a complex refractory residue that might account for the dark, non-icy material on the surface (see review in Cruikshank et al., 2015).
Finally, the last hypothesis identifies the dark red material covering the low reflectance terrains of Pluto as haze particles deposited onto the surface. A layered haze extending up to about 350 km above the surface has been discovered by the New Horizons spacecraft during its flyby (e.g. Gladstone et al., 2016;Grundy et al., 2018;Young et al., 2018). Haze particles originate from the dissociation and ionisation of Pluto's upper atmosphere gases (CH 4 , N 2 , and CO) due to radiolysis and photolysis processes. Different energetic sources are possible drivers for the haze production (e.g. Cheng et al., 2017;Grundy et al., 2018). First, photons originating from the Sun are an important energy source. Ly-photons are energetic enough to break C-H bonds in CH 4 and their fluxes are relatively high. Although more energetic photons are present in smaller quantities than Ly-photons, they also contribute to photolysis and are able to break the stronger triple bonds in N 2 and CO gases. In addition to photons, protons and electrons originating from the Sun are also expected to trigger gases dissociation in Pluto's upper atmosphere. Finally, cosmic rays also play a role despite their low fluxes at Pluto's distance, as they deliver highly energetic particles.
The formation of haze particles starts in the upper layers of Pluto's atmosphere, where dissociation products (ions, radicals) accrete and build larger and more complex molecules. Other gases like HCN, C 2 H 2 , or H 2 O can stick to the haze particles, either being adsorbed, chemically bonded, or directly condensed around the haze particle (e.g. Luspay-Kuti et al., 2017;Gao et al., 2017;Wong et al., 2017). When growing, the haze particles also become heavier and start falling down. The temperature rises as they reach lower altitudes, initiating sublimation of the most volatile constituents of the haze particles.
The analysis of combined MVIC and LEISA over both Cthulhu and Lowell regions by Protopapa et al. (2020) showed that a single colouring agent could account for the diversity of Pluto's surface colours. The different colours are likely to be caused by changes in abundance and particle size of a single material, quite similar to Khare's Titan tholins. Although this does not discard the other scenarios, it supports the hypothesis of an atmospheric origin for Pluto's dark material.
This paper further investigates this hypothesis, by analysing whether analogues of the haze particles forming in Pluto's atmosphere could reproduce the reddest colours observed on the surface. The spectral reflectance properties of the analogues (so-called tholins, as introduced in Sagan and Khare, 1979) were measured through spectro-gonioradiometry, and were then fitted to the spectra collected by the Ralph instrument of the New Horizons spacecraft during its flyby of Pluto.
The New Horizons spectrophotometric observations are first presented in Section 2, while the experimental setups used to synthesise the analogues of Pluto's haze particles and measure their reflectance spectra are described in Section 3. The numerical models used to simulate the surface reflectance are then detailed in Section 4. The experimental results are provided in Section 5, where the optical properties of the Pluto tholins are analysed. Section 6 presents the comparison with respect to New Horizons data. Finally, results are discussed in Section 7 and the main concluding points are summarised in Section 8.

New Horizons observations
Pluto's flyby performed by the New Horizons spacecraft occurred in July 2015, the point of closest approach having been reached on July, 14th after about 10 years of travel in space (Stern et al., 2015). The composition of Pluto's surface was investigated through observations collected by the Ralph instrument (Reuter et al., 2008). It combines a visible and near-infrared camera (MVIC) and an infrared imaging spectrometer (LEISA).

Measurements
MVIC possesses four spectral channels, the first three being referred to as the BLUE filter (covering a spectral range extending from 0.40 μm to 0.55 μm), the RED filter (from 0.54 μm to 0.70 μm), and the NIR (near-infrared) filter (from 0.78 μm to 0.98 μm), respectively. Additionally, one channel covering wavelengths from 0.86 μm to 0.98 μm is dedicated to the detection of the narrow CH 4 band present in this specific spectral range.
The LEISA instrument covers a spectral range going from 1.25 μm to 2.50 μm. Most of the chemical components expected to be present on Pluto's surface show absorption features in this spectral range, thus making LEISA data able to verify their identification. The spectral resolution of the LEISA instrument is ∕ = 240 between 1.25 μm and 2.50 μm. A higher spectral resolution ( ∕ = 560) was supposed to be available from 2.1 μm to 2.25 μm, but this higher resolution segment turned out to be difficult to calibrate properly. Therefore, this work focuses on the lower spectral resolution segment only.
The maps generated by the LEISA instrument comprise a spectral dimension and a spatial one. The near-infrared spectrum from 1.25 μm to 2.5 μm is obtained along the in-track direction. The reflectance of different parts of the surface is theoretically measured at the same wavelength in the cross-track direction. By scanning the image in the in-track direction, and thus at different wavelengths, the full spectrum of a given region of Pluto can be constructed. However, due to slight instrumental effects, the measurements are not exactly conducted at the same wavelength in the cross-track direction. A so-called ''spectral smile'' of a few spectels is present over the 256 pixels array along the cross-track direction. During the scanning of the target by LEISA, the spacecraft kept rotating around its -axis so that Pluto's surface is scanned in the transmission band pass direction (in-track direction). It also tries to compensate any shift of the scan in the perpendicular direction with brief thruster bursts. These zig-zag scans have to be corrected during data projection.

Data pre-processing
Both photometric and spectral calibrations of the Ralph instrument were conducted before the launch of the spacecraft and completed with in-flight calibrations with respect to stars. An additional ''flatfield'' calibration was performed, based on a rapid scan of Pluto's surface in the cross-track direction (as opposed to the in-track direction which is the nominal orientation of the instrument when collecting observations). Each pixel of the 256 pixels array thus covers almost the same part of the surface. The overall calibration of LEISA data is based on a combination of those ground-based and in-flight stellar calibrations, spectral smile corrections, and flat fielding. Nonetheless, some calibration issues remain in New Horizons data. The LEISA reflectance values are overestimated by an amount of 26% and should therefore be scaled down accordingly (Protopapa et al., 2020).
Several scans of Pluto's surface by the LEISA instrument were conducted, and data cubes were reconstructed with the Integrated Software for Imagers and Spectrometers (ISIS) software from the United States Geological Survey (USGS). Two of the close-encounter data cubes were calibrated, and projected onto LORRI Pluto-sized spherical data map (orthographic projection at same observation geometry). Sharp geological features identifiable in both LORRI and LEISA data were finally used for slight mismatch corrections in LEISA data projection, with respect to the higher resolution LORRI base map. Together, the two calibrated LEISA data cubes cover the whole visible disc of the planet.
In this work, we used the same calibrated data cube covering Cthulhu region as in Schmitt et al. (2017). This data cube does not consist of raw data but follows from a Principal Component Analysis (PCA), described in more details in Schmitt et al. (2017). This postprocessing step aims at identifying the main surface components, as well as mapping their distribution, and removing most of the instrumental effects and noise. The result of the PCA was that the first nine Principal Components account for almost 95% of the variance, demonstrating the good quality of LEISA data. Reduced-noise data cubes were reconstructed from a PCA inversion, by assuming that the first Principal Components account for the major part of the variance. Using these reduced-noise data cubes instead of the raw ones allows for a sharper analysis of New Horizons spectra collected over the Cthulhu region.
The LEISA data from the reconstructed cube covering Cthulhu region were scaled down by 26% to correct the calibration issues previously mentioned. A precise and reliable value for the uncertainty of this calibration factor is hard to determine. We however used 7% as a conservative estimate for MVIC data error bars. This value is based on the comparison of two independent calibration methods for MVIC data in Howett et al. (2017). Discrepancies in the correction factor are limited to 3%-7%, such that 7% seemed a reasonable upper limit for the uncertainty value. Finally, MVIC data were corrected to LEISA geometry, following the method presented in Protopapa et al. (2020). New Horizons reflectance spectra. The light and dark blue curves correspond to New Horizons LEISA spectra collected over H 2 O-rich and H 2 O-poor eastern regions of Cthulhu, respectively. They are retrieved from the reconstructed data cube described in Section 2.2. MVIC data are corrected to LEISA geometry (see Protopapa et al., 2020). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Average coordinates and viewing geometry of the two regions of interest for this study (H 2 O-poor and H 2 O-rich regions of Cthulhu), over which the New Horizons spectra displayed in Fig. 1 were retrieved and averaged.

Cthulhu spectra
To identify the red dark material covering Cthulhu, New Horizons spectra collected over this specific region were extracted from the data cube described above. The main focus of this study is on the core of the Cthulhu region, where the dark organic material is supposed to be the most concentrated. A spatial average was then computed over several pixels which were identified as the end-member of the red material in Schmitt et al. (2017). Using the averaged spectrum as a reference for comparison with laboratory reflectance measurements allows to further reduce the pixel-to-pixel noise of the New Horizons spectrum.
The MVIC data and LEISA reflectance spectrum of the red material's end-member are reproduced in Fig. 1. In addition to the dark material end-member spectrum, we also used an average spectrum collected in the eastern periphery of Cthulhu. This spectrum exhibits strong absorption bands at 1.5 μm and 2.0 μm, demonstrating the presence of water ice mixed with the dark material. Evidences of a gradient in H 2 O concentration from the centre to the periphery of Cthulhu were already brought in Schmitt et al. (2017).
Because the reflectance level is highly dependent on the illumination and observation geometry under which a spectrum is collected, the average geometries corresponding to our two spectra of interest were retrieved. They are reported in Table 1, along with the average coordinates of the two regions over which these spectra were obtained.

Experimental setup
The way Pluto tholins are synthesised is first described in Section 3.1, before presenting the experimental setup used to conduct the laboratory reflectance measurements in Section 3.2.

Tholins synthesis
Pluto tholins were synthesised in LATMOS (French acronym for Laboratoire Atmosphères, Milieux, Observations Spatiales), with the PAMPRE experimental setup (Production d'Aérosols en Microgravité par Plasma REactif, for more details see e.g. Szopa et al., 2006;Alcouffe et al., 2009;Jovanović et al., 2020). It was first designed to synthesise laboratory analogues for Titan's haze particles, from an initial gas mixture of N 2 and CH 4 . For this study, it was used to produce analogues of Pluto's haze particles, by adding CO to the N 2 ∶CH 4 gas mixture (Lellouch et al., 2017;Young et al., 2018). PAMPRE simulates the dissociation of gases molecules as it occurs in Titan's or Pluto's high atmospheres, due to UV light or energetic particles. In the PAMPRE setup, the molecules are dissociated in a cold plasma discharge between two electrodes. Dissociation products evolve into more complex products. As they grow, they become heavier and are finally recovered as a dusty material made of small spherical particles. The synthesis is conducted at low pressure (0.9 ± 0.1 mbar) and room temperature.
To synthesise Pluto tholins, two initial gas compositions were used: N 2 ∶CH 4 =99%:1% and N 2 ∶CH 4 =95%:5%, with 500 ppm of CO being added to each of them. They are thought to be representative of Pluto's atmosphere composition at altitudes of 400 and 650 km, respectively (Lellouch et al., 2017;Young et al., 2018). Their chemical composition was studied in Jovanović et al. (2020), and their optical properties in Jovanović et al. (2021).
In the following, for the sake of conciseness, these two types of Pluto tholins are referred to as 5%CH 4 and 1%CH 4 tholins, since the relative amount of CH 4 drives the reflectance differences between them (see discussion in Section 6.1.2). Small quantities of these Pluto tholins were also mixed with pyrrhotite, an opaque mineral used here as a darkening agent to model the effect of contamination by dark interplanetary dust (see discussion in Section 7).

Laboratory reflectance measurements
The spectrophotometric properties of the tholins samples were investigated by measuring their reflectance at various geometries with a spectro-gonio-radiometer developed at IPAG: the Spectro-photometer with cHanging Angles for Detection Of Weak Signals (SHADOWS) . This instrument operates in the 0.4-4.8 μm spectral range and thus covers the full spectral range of the New Horizons observations both in the visible and in the near-infrared (from 0.4 μm to 2.5 μm), with an absolute photometric accuracy better than 1%. It therefore allows for a direct comparison between laboratory measurements and spectra collected by the New Horizons spacecraft.
The spectro-gonio-radiometer can perform reflectance measurements under various illumination and observation geometries. The geometry under which the spectrophotometric measurements are conducted significantly affects both the reflectance level and the spectral contrast (Hapke, 2012b). To compare the spectra measured in the laboratory on synthesised samples with those collected by New Horizons, the illumination and observation geometries must therefore be as similar as possible. This is necessary to remove any variation of the reflectance induced by a geometry difference and not by the reflectance properties of the material itself.
Nonetheless, first measuring the reflectance spectra of the tholins samples under a larger range of illumination and observation geometries allows to characterise the photometric function of the samples, as well as the evolution of their absorption features as a function of the geometry. This is crucial to retrieve tholins' optical properties from their reflectance measurements for modelling purposes (see Section 5.5). The reflectance measurements were conducted under incidence angles of 0 • , 30 • , and 60 • , and emergence angles ranging from −70 • to 70 • with a step of 10 • . Some additional measurements were carried out under vacuum and with progressive heating of the samples (up to 100 • C), to verify that the observed absorption features were not caused by adsorbed water.

Numerical models
Complementing the laboratory reflectance measurements, an analytical model was used to simulate the reflectance spectra of Cthulhu's surface. Section 4.1 first discusses the reflectance model, while Section 4.2 describes the modelling of the surface itself.
4.1. Reflectance model 4.1.1. Nominal Hapke model In this work, the Hapke model was selected, as one of the most complete theoretical model currently available to simulate surfaces' reflectance (e.g. Hapke, 1986Hapke, , 1993Hapke, , 2002Hapke, , 2008Hapke, , 2012a. The complete formulation of the Hapke model is given as follows: with 0 = cos( ) and = cos( ). The reflectance is given as a function of the illumination and observation geometry angles , and (incidence, emergence, and phase angles respectively), and of various parameters describing the surface properties. The angles and refer to the incidence and emergence angles corrected for roughness. The porosity and macroscopic roughness of the surface are both taken into account (parameters and̄, respectively). The shadowing function describes how shadows are cast on the surface particles, and the resulting change in the total reflectance. It accounts for the effect of the macroscopic roughness̄at a given illumination and observation geometry. The brightness growth observed at very small phase angles is described by the two opposition effect terms: SH and CB , representing the shadow hiding and coherent backscattering opposition effects, respectively. Multiple scattering is accounted for in the function. Finally, the phase function ( , ) describes the dependency of the photometric level to the phase angle .

Simplified reflectance model
Additional assumptions were made to simplify the complete formulation of the Hapke model given in Eq. (1). Over the Cthulhu region, New Horizons spectra were collected at phase angles close to 22 • (Table 1). The contribution of the coherent backscattering opposition effect term in the reflectance model can thus be safely neglected. Indeed, this effect generally only manifests itself for phase angles smaller than 2 • (Hapke, 2002).
The multiple scattering function is still included in our simplified version of Hapke model (analytical expressions for can be found in Hapke, 2012b). The reflectance measurements performed on Pluto tholins indeed showed that they are extremely bright in the nearinfrared. Multiple scattering cannot be neglected in that case since the absorption might not be efficient enough to ensure that the incident light is fully absorbed after encountering a single particle.
Finally, a macroscopic roughness of 20 • was assumed for the shadowing function ( , , ,̄), which is the value derived from Hapke model's fit of the phase curve obtained during 2018 Pluto's occultation (Verbiscer et al., 2019).

Surface constituents
The low reflectance terrains of Pluto are not only covered by a dark red organic material. As already discussed in Section 2.3, water was identified towards the periphery of Cthulhu, where the dark material is probably mixed within H 2 O ice. CH 4 was also detected on the topographic heights of Cthulhu and could correspond to seasonal deposits only subsisting at high altitudes. CH 4 ice could also be present in the bedrock of Pluto, along with H 2 O which is a potential bedrock constituent . Additionally, tentative evidences for hydrocarbon-rich ices in Cthulhu region were brought in Cook et al. (2019). Hydrocarbon-rich species were also mentioned as photochemical products participating in the formation of haze particles in Pluto's atmosphere in Grundy et al. (2018).
The model of Pluto's surface proposed in this work is thus a mixture of the above-mentioned icy species, namely H 2 O, CH 4 , as well as CH 3 OH and C 2 H 6 which were selected as representative for the potential hydrocarbon species present in Cthulhu. Pluto tholins were added, as analogues for the dark-red organic material.

Intimate mixing between refractory species
A simple mixing model representing the surface as a spatial mixing between two separate units was first considered (represented in Fig. 2). The first unit contains the volatile species (here CH 4 ) while the second one is composed of the refractories (tholins, H 2 O, CH 3 OH, and C 2 H 6 ).
The total reflectance of the surface is calculated as the sum of the reflectance of the two units, weighted by their relative spatial ratio (left as an unknown in our model). Within the refractory unit, it was assumed that the different species are mixed at the individual particle scale. The reflectance of such an intimate, homogeneous mixture between various chemical compounds is modelled from an averaged single scattering albedo and phase function. These averaged values are computed from the individual properties of the different compounds involved in the mixing (relative mass ratio , density , particle mean diameter and single scattering albedo ), as follows:

Ices condensed on tholins particles
A second mixing model was also considered, which is thought to be more representative of the way haze particles form in Pluto's atmosphere. As haze particles grow and fall down to the surface, ices can condense around them (e.g. Luspay-Kuti et al., 2017;Gao et al., 2017;Wong et al., 2017). The particles then sediment and settle onto Pluto's surface, forming an intimate mixture of tholins surrounded by various ices.
With this surface model, tholins are thus not directly mixed with icy components, but rather trapped into icy particles which are intimately mixed with each other. Tholins act as a colouring pigment in an icy matrix. To simulate the reflectance of ices condensed around tholins particles, Maxwell-Garnett effective medium model was applied. The use of this model is justified by the fact that the tholins particles are smaller than the wavelength. The effective dielectric constant ef f can be deduced from those of the matrix and of the inclusions, and is given by where is the fraction of the embedded material among the matrix, while 1 and 2 refer to the dielectric constant of the matrix and of the embedded material, respectively. The effective optical constants of such a medium can be directly retrieved as the real and imaginary parts of the dielectric constant, as follows: The effective optical constants were thus computed for tholins surrounded by various ices, and an intimate mixing was assumed between the different condensed ice particles. Still, a spatial mixing was modelled between one spatial unit covered by intimately mixed H 2 O, C 2 H 6 and CH 3 OH icy particles condensed around core tholins, and a second unit composed of CH 4 -ice to account for the possible presence of CH 4 -ice in the bedrock of Pluto. The spatial mixing ratio between the two units is again an unknown of the model. A schematic representation of this second mixing model is provided in Fig. 3.
An additional difficulty arises from the use of the Maxwell-Garnett model. While the simplified Hapke model only requires the single scattering albedo of the Pluto tholins tholins , the Maxwell-Garnett model needs the associated optical constants tholins and tholins to compute the effective optical constants (see Eqs. (10)- (12)). The single scattering albedo of the effective medium is then computed from ef f and ef f and plugged in the Hapke model to compute the reflectance. However, assuming tholins is known, there is no direct way to uniquely determine and from the single scattering albedo, as clearly shown in Eqs. (7) and (8).
A commonly used trick is to assume constant and find the wavelength-dependent values of that fit the single scattering albedo. However, a constant value for tholins over our spectral range of interest seems rather unrealistic considering the strong absorption in the visible, responsible for the observed red slope. An alternative option is thus to use a set of values available in the literature which reproduces the photometric level of the tholins reasonably well. Using those values, the values are computed to match the single scattering albedo, relying on the fact that most of the spectral information is carried by . This second method is expected to provide more realistic estimates of tholins optical constants and was therefore applied in this work, using the values of provided in Khare et al. (1984) (see Section 5.4 for comparison of different sets of tholins optical constants). More details about the numerical determination of the optical constants and from are given in Appendix A.2.

Fitting the model to the New Horizons observations
We tried to reproduce New Horizons spectra using the models presented in Sections 4.1 and 4.2. The surface was modelled as a mixture between H 2 O ice, hydrocarbon-rich ices (e.g. C 2 H 6 , CH 3 OH), CH 4 ice and tholins, as analogues for the dark material covering Cthulhu region.  CH 3 OH-ice: Trotta (1996) and Cook et al. (2019). Most of these data are available online in the GhoSST database of the SSHADE infrastructure (www.sshade.eu).
For simplification purposes, the single scattering phase function of Cthulhu's surface was assumed to be isotropic and independent of wavelength (i.e. ( , ) = 1). This was mainly motivated by the fact that the phase functions of the surface constituents are unknown.
When computing the reflectance of the modelled surface, many parameters remained unknown, among which the relative abundance and mean particle diameter of each surface constituent. The macroscopic porosity parameter was also left as a free parameter. An optimisation algorithm was thus applied in order to find the best fit with respect to the New Horizons spectra.
The final objective was to assess the relevance of Pluto tholins as analogues of the dark material covering Cthulhu region. It was therefore of primary importance to not only match the near-infrared part of the spectrum but also MVIC data, collected in the visible. The darkest regions of Pluto indeed present a very strong red slope in the visible which is due to the dark material only, since ices have very high photometric level in both the visible and the near-infrared. The fit was thus determined by computing the 2 with respect to both MVIC and LEISA data. Because the MVIC instrument only sampled three different wavelengths and therefore had a much poorer spectral resolution than the near-infrared spectrometer, MVIC data were weighted to give them as much weight as the whole LEISA spectrum, strongly constraining Cthulhu spectra in the visible.
Different optimisation algorithms were tested and tuned to find the best performing one, both in terms of convergence and computational effort. A simple genetic algorithm was eventually selected as offering the best trade-off between these two criteria for our particular problem (see Appendix B for further details on the optimisation process).
It must be stressed that the global optimum was neither targeted nor achieved with the optimisation process described above. The number of parameters to estimate and the search space were both too large, and critical information was missing to precisely reproduce the reflectance of Cthulhu's surface. However, the main interest of this approach was not to estimate the exact amount of tholins and/or icy species present in the Cthulhu region. It was rather to determine whether modelling the reflectance of tholins mixed with different ices could help reproducing the spectra collected by the New Horizons spacecraft. This also made our simplifying assumptions acceptable here (e.g. isotropic phase function).

Experimental results
The size and shape of Pluto tholins were first characterised from Scanning Electron Microscopy (SEM) observations. The laboratory reflectance measurements were then analysed, to better understand the reflectance properties of our tholins and thus make the comparison with New Horizons spectrophotometric data possible.

SEM observations
The SEM analysis was performed on the two types of Pluto tholins (5%CH 4 and 1%CH 4 , respectively), with the JEOL JSM-840 A SEM at University Paris VI and the ZEISS Ultra55 FEG-SEM at CMTC INP Grenoble operating at 10 kV. It revealed tholins are very regular spherically-shaped particles, which sometimes agglutinate to form larger aggregates. The difference in initial gas mixture composition leads to some disparities in the physical properties of the synthesised particles. 5%CH 4 tholins overall show more regular shapes than 1%CH 4 ones. They appear as almost spherical, while the 1%CH 4 tholins particles present more irregularities, still being spherically-shaped (see Fig. 4).
From SEM observations, the grain size distribution of Pluto tholins was estimated by measuring the apparent radius of tholins particles, averaged over a few hundred particles for each of the two tholins compositions (see Fig. 5). The mean radius of both 5%CH 4 and 1%CH 4 tholins particles was estimated equal to 210 nm although the median is lower, especially for 1%CH 4 tholins (around 200 nm). The grain size distribution of 5%CH 4 tholins is slightly narrower (standard deviation = 35 nm) than that of the 1%CH 4 tholins ( = 55 nm), an observation already made for Titan tholins (Hadamcik et al., 2009). The more irregular shapes of 1%CH 4 tholins might introduce some additional uncertainty in the mean diameter estimation, but cannot fully account for the much wider dispersion of the mean radius values. The difference in the dispersion of the grain size distribution is clearly observable even without performing any statistical analysis from the SEM images (see Fig. 4).

Reflectance spectra as a function of the geometry
The influence of the illumination and emergence geometry is clearly visible from our laboratory reflectance measurements. Fig. 6 displays the reflectance spectra of 5%CH 4 tholins in the 0.4-2.5 μm range. These spectra were obtained under a constant incidence angle set to 0 • , but with varying emergence angles (ranging from 10 • to 70 • ), and thus with changing phase angles. The photometric level strongly reduces with increasing phase angle. However, the absorption features are only slightly affected by the geometry variations, becoming shallower when the phase angle increases. Fig. 7 also illustrates the evolution of the tholins reflectance with varying incidence and emergence angles. For each of the three incidence values (0 • , 30 • and 60 • , respectively), the reflectance is displayed as a function of the emergence angle, for four different wavelengths (all comprised in the New Horizons spectral range). The tholins tend to reflect the incident light backwards, which is in agreement with the fact that the reflectance decreases with increasing phase angle. Such a backscattering could possibly be due to the very small size of the tholins particles (comparable to the wavelength). From both Figs. 6 and In both samples, tholins are submicrometer-sized spherical particles. 1%CH 4 tholins however present more irregularities and show more diverse shapes compared to 5%CH 4 tholins. 1%CH 4 tholins also seem to present a wider range of sizes.

Fig. 5.
Grain size distribution derived from SEM observations. Both distributions have an average radius of about 210 nm, although the median is slightly lower for 1%CH 4 tholins (around 200 nm). The size distribution of 1%CH 4 tholins is broader, with a standard deviation of = 55 nm, compared to that of 5%CH 4 tholins ( = 35 nm). 7, it is clear that the observation geometry effects need to be accounted for, and that the comparison between laboratory measurements and New Horizons data is only meaningful if conducted at similar incidence, emergence and phase angles.
Figs. 6 and 7 only display results obtained with the 5%CH 4 tholins but identical reflectance measurements were conducted with the 1%CH 4 sample. There is no noticeable difference between the two types of Pluto tholins when looking at the photometric behaviour (1%CH 4 tholins are also characterised by a strong backscattering behaviour). However, the composition of the initial gas mixture of the tholins synthesis has a strong influence on the shape and absolute level of the reflectance spectra. Both 5%CH 4 and 1%CH 4 tholins' spectra are provided in Fig. 9 (geometry: = 0 • , = 30 • ). Reducing the relative amount of CH 4 in the synthesis gas darkens the resulting tholins, in agreement with previous findings for Titan tholins . 1%CH 4 tholins show a slightly lower reflectance level in the near-infrared. They also exhibit a more progressive increase of reflectance in the visible, but a steeper slope in the very nearinfrared (1.0-1.3 μm). These differences between the two types of tholins become of particular interest when comparing the laboratory reflectance spectra with New Horizons data (see Section 6).

Photometric and spectral effects of mixing tholins with a darkening agent
The deposition and mixture of opaque interplanetary dust grains with the dark red organic material could darken and change the spectral contrast of Pluto's surface with time. It is thus crucial to quantify the effect of mixing Pluto tholins with a darkening agent on both the photometric level and spectral constrast. To this end, the reflectance measurements of Pluto tholins mixed with pyrrhotite (introduced in Section 3.1) are provided in Fig. 8.

Comparison with optical constants of tholins from literature
We investigated whether optical constants derived for other tholins samples and already available in the literature could be used to reproduce the reflectance of the Pluto tholins analysed in this work. Various types of tholins were previously synthesised and analysed in the literature, mostly to characterise Titan aerosols (Brasse et al., 2015). Multiple techniques exist for the tholins synthesis (e.g., Coll et al., 2013), each of them associated with very different experimental conditions (various temperature and pressure conditions, energy source, gas flow rate, irradiation duration, etc.).
The final composition of the tholins particles, and therefore their optical constants, depend on the technique employed to synthesise them. Additionally, the method used to measure the optical constants also influences their values (Brasse et al., 2015). For these reasons, optical constants of tholins available in the literature are very different from one sample and from one measurement technique to another.
An attempt was made to compare the reflectance measurements performed on our N 2 ∶CH 4 ∶CO tholins with spectra calculated from various N 2 ∶CH 4 Titan tholins' optical constants available in the literature (Brasse et al., 2015). Results of this comparison are displayed in Fig. 9. It is especially interesting to note that very diverse slopes are obtained in the visible, even between tholins synthesised from initial gas mixtures with similar composition. Unfortunately, none of the available optical constants sets led to a proper fit of our laboratory spectra.
Optical constants presented in Mahjoub et al. (2012) and Sciamma-O'Brien et al. (2012) were determined from tholins synthesised under similar conditions and using the same experimental setup as the Pluto tholins analysed in this study. But, even these optical constants measured on rather similar samples could not reproduce the reflectance spectra obtained from our laboratory measurements. The simulated reflectance spectra obtained with the optical constants of their 5%CH 4 tholins show a steeper red slope in the visible, which is also shifted towards lower wavelengths compared to our laboratory measurements. On the other hand, the slope of their 1% CH 4 tholins is much weaker. Moreover, the reflectance level of our Pluto tholins is much higher in the near-infrared. It must however be stressed out that the synthesis process was significantly different between our Pluto tholins (powder material synthesised while in suspension in the plasma discharge, see Section 3.1) and the Titan tholins analysed in the above-mentioned studies (organic films deposited on various substrates). On the contrary, optical constants presented in Khare et al. (1984) reproduce the 5%CH 4 tholins' visible spectral slope and very nearinfrared photometric level reasonably well (see Fig. 9). These optical constants are the most commonly used, but they have a poor spectral resolution. Despite the relatively good fit in the 0.4-1.3 μm range, the four absorption bands observed above 1.3 μm in our laboratory spectra do not appear in the simulated spectra obtained with the optical constants found in Khare et al. (1984).

Retrieving single scattering albedo from reflectance measurements
The previous analysis of existing sets of optical constants clearly showed that it is not possible to reproduce our laboratory reflectance measurements from optical constants found in the literature. Consequently, the single scattering albedo of the Pluto tholins needs to be determined in a different way. It is possible to isolate and estimate it from the many reflectance spectra collected under various geometries. Indeed, the phase function ( , ) and the single scattering albedo ( ) are the only unknowns in the simplified Hapke model formulation given by Eq. (1).
While the single scattering albedo depends on the wavelength only, the phase function depends on both the wavelength and the phase angle . The wavelength-dependency of the phase function is however relatively weak (Hapke, 1993(Hapke, , 2012a and was thus neglected in this analysis. The single scattering albedo therefore accounts for the full spectral information, while the phase function ( ) is a function of the illumination and observation geometries only.
The laboratory reflectance measurements collected under different illumination and observation conditions, and over a wide spectral range, were all combined together. A non-linear least-squares inversion of the simplified Hapke model was then applied to retrieve the functions ( ) and ( ). There was no more correlation between the single scattering albedo and the phase function, and therefore they could be estimated simultaneously. Because the reflectance model is highly nonlinear, the inversion required an iterative process, converging towards an unique solution by successive differential corrections from a given initial guess. Some convergence issues were first encountered because of the intrinsic constraints of the single scattering albedo values (ranging from 0 to 1). A numerical singularity indeed exists when tends towards its maximum value, and deteriorates the convergence of the least-squares algorithm. An extended grid search simulation was thus run to identify a good initial guess to initiate the least-square algorithm. The results of the Hapke's model inversion are presented in Figs. 10 and 11. The consistency between the grid-search and least-squares solutions for the single scattering albedo is a good indicator of the convergence of the inversion. The solutions obtained for the phase function for both the 5%CH 4 and 1%CH 4 tholins are also insightful. While the phase functions generated by the grid-search analysis are different for the two tholins, the converged least-squares solutions are conversely extremely similar (see Fig. 11). Such significant differences between the phase functions of the two types of tholins obtained from the grid-search analysis solutions are not expected, clearly showing the limitations of the initial guess provided by the grid-search. The results on the phase function also highlight the good convergence of the hybrid inversion technique and demonstrate the improvement brought by the least-squares algorithm compared to the coarse grid-search solution.
The single scattering albedo of Pluto tholins exhibits a very steep red slope in the visible, which seems consistent with the New Horizons observations collected over the dark-red terrains of Pluto. In the near-infrared, is very close to 1, in agreement with our laboratory reflectance spectra which show tholins as extremely bright in the nearinfrared. It also explains the difficulties encountered by the simple least-squares algorithm to converge to an unique solution so close to the numerical singularity if not combined with a prior grid search analysis.
Furthermore, most of the spectral information seems to be captured by the inverted single scattering albedo. This also proves that the inversion was successful, since was the only parameter assumed to be wavelength-dependent. The four absorption bands identified in our laboratory spectra are indeed present in the inverted single scattering albedo. Finally, the phase function is globally decreasing with increasing phase angle. This finding is in agreement with the strong backscattering observed in Fig. 7. For both the 1%CH 4 and 5%CH 4 tholins, the inverted single scattering albedo values are reported in Table A.1 (as a function of the wavelength) and the phase functions can be found in Table A.2 (see Appendix A.3).

Comparison with MVIC/LEISA observations
Using Pluto tholins' reflectance properties as described in Section 5, their relevance as analogues for the dark red material covering Cthulhu region was assessed by comparing their reflectance spectra with the M. Fayolle et al. New Horizons data collected by the Ralph instrument when flying over Cthulhu.

Reflectance spectra
Before trying to reproduce the New Horizons data using numerical models exposed in Section 4, Cthulhu spectra were directly compared with laboratory reflectance measurements of Pluto tholins. The comparison was performed at identical illumination and observation geometries, to prevent any bias in the photometric level caused by a difference in the measurement geometries. The laboratory spectra were convolved to match the spectral resolution of the New Horizons observations, which in turn were corrected to account for photometric calibration issues. Fig. 12 shows the superposition of Cthulhu reflectance spectra with those of the pure tholins samples, under the same geometry. Despite a good agreement at first order for the 1%CH 4 tholins, there are striking differences between tholins' and New Horizons' spectra, both spectrally and photometrically. The four absorption bands in tholins laboratory spectra (at 1.5, 1.75, 2.0 and 2.3 μm) are completely absent in the New Horizons spectra, and the same holds for the spectral features observed in New Horizons data which do not appear in the tholins spectra.
In addition, Fig. 12 clearly highlights the influence of the tholins composition. In the near-infrared, the reflectance of 1%CH 4 tholins is much closer to New Horizons spectra in terms of photometry (being about 20% lower than that of 5%CH 4 tholins on average). The red slope displayed by the 1%CH 4 tholins in the visible is also more consistent with MVIC data.

Absorption bands in laboratory reflectance spectra
The reflectance spectra of the tholins show several absorption features in the near-infrared, as displayed in Fig. 13. As mentioned above, four bands clearly appear in the laboratory measurements at about 1.5 μm, 1.75 μm, 2.0 μm and 2.3 μm.
Additional reflectance measurements were conducted on COdepleted N 2 ∶CH 4 tholins to quantify how the presence of a small amount of CO in the initial gas mixture affects the tholins spectrophotometric properties. The result of this comparative analysis is displayed in Fig. 14, but no significant spectral difference and/or band shift can be observed between the tholins with and without CO. The small difference in photometric level can be at least partly attributed to slightly different properties of the samples (grain sizes, density, compaction, etc.). These results indicate that the presence of CO is not responsible for any additional spectral signature in the tholins reflectance spectra. The absorption bands observed in tholins spectra are thus more likely due to bonds between C, H and N compounds rather than to bonds involving oxygen.
The band around 1.75 μm appears to be deeper and wider for 5%CH 4 than for 1%CH 4 tholins. It can be assigned to the first overtone of the [C-H] bond (Workman, 1996). This is in agreement with the fact that 5%CH 4 tholins are expected to contain more [C-H] bonds, explaining the stronger absorption of this feature for 5%CH 4 compared to 1%CH 4 tholins.
The absorption feature around 2.25-2.3 μm is thought to be related to the first overtone of the stretching mode of [C≡N] bonds, as proposed in Cruikshank et al. (1991). Another possible interpretation for this band was proposed by Cloutis (1989), assigning it to combination of aromatic carbon stretch with [C-H] stretching modes, and combinations of [C-H] stretching and bending modes. However, the shape and depth of this band does not seem to strongly depend on the N 2 /CH 4 ratio of the initial gas mixture from which the tholins are synthesised, on the contrary to what has been observed for the 1.75 μm band. The fact that this 1.75 μm absorption feature conserves its shape, width and depth in the reflectance spectra of the two types of tholins thus appears more consistent with its assignment to [C≡N]

bonds rather than to [C=C] and [C-H] combinations.
The two remaining bands, located around 1.5 μm and 2.0 μm respectively, could at prima facie correspond to H 2 O absorption features. However, as mentioned in Section 3.2, additional laboratory reflectance measurements were conducted under vacuum and after heating to verify that adsorbed water was not responsible for any absorption band. The 1.5 μm band is thus more likely related to the first stretching overtone of [N-H] bonds, while the 2.0 μm band can be due to some combinations of [N-H] bending and stretching modes (Workman, 1996). Nevertheless, it must be noted that the presence of those amine functions could be partially related to the synthesis experimental setup, as the formation of nitrogen-bearing compounds is enhanced by wall effects.

The assignment of some tholins' near-infrared absorption bands to [C-H] and [N-H] functions is consistent with previous findings
showing that these functions contribute to N 2 ∶CH 4 tholins' mid-infrared spectra (Quirico et al., 2008;Gautier et al., 2012). Comparative reflectance measurements of tholins with and without CO∶5%CH 4 (or 1%CH 4 ) tholins without CO compared to 5%CH 4 (or 1%CH 4 ) tholins with 500 ppm CO (for CO−depleted tholins, measurements were only conducted from = 0.68 μm). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Absorption bands in New Horizons spectra
The New Horizons spectra show two main absorption features (see Figs. 1 or 12). The first one is located around 2.1 μm. When getting closer to Cthulhu periphery where water ice exists in larger quantities, the H 2 O band is contributing to this absorption feature. However, this band still appears in the centre of Cthulhu (i.e. in the reddest terrains of Pluto), where very little water ice is expected (as other H 2 O absorption bands are not present there). The interpretation of this band is still not firmly resolved (Schmitt et al., 2017).
In addition, another absorption feature is observed around 2.3 μm in the Cthulhu region. Its presence is however thought to be related not to the dark organic material itself, but rather to some hydrocarbon ices with which the red material is expected to be mixed (Cook et al., 2019). Oddly, the shape and depth of this band remain unchanged from Cthulhu's centre to its periphery. The presence of more water ice on the borders of Cthulhu would suggest that the red material is more diluted in water ice there, and more concentrated in Cthulhu's centred, reddest parts. The fact that the 2.3 μm band seems independent of the relative amounts of dark material and water ice favours the assignment of this band to hydrocarbon ices, as suggested in Cook et al. (2019). Furthermore, both Cthulhu spectra exhibit two additional weaker bands at about 1.3-1.35 μm and 1.8-1.85 μm, which might also be linked to these hydrocarbon ices. Fig. 15 shows near-infrared spectra of possible hydrocarbon-ices candidates C 2 H 6 , C 3 H 8 and CH 3 OH, superimposed upon New Horizons Cthulhu spectra.

Optimised fitting between New Horizons data and simulated reflectance
The optimisation process described in Section 4.3 was applied to investigate whether the New Horizons spectra could be reproduced by combining the laboratory reflectance measurements of Pluto tholins (Section 5.2) with the numerical models developed in Sections 4.1 and 4.2.
To account for the diversity of the New Horizons spectra measured over the dark terrains of Pluto, we tried to reproduce the reflectance behaviour of both Cthulhu's H 2 O-rich region and H 2 O-poor regions, as discussed in Section 2.3. The latter corresponds to the reddest terrains of Pluto and is thus expected to contain the largest relative amounts of the dark Plutonian material.

Intimate mixing between refractory species
The fitting process was first applied to the simple intimate mixing surface model presented in Section 4.2.2. For each species, both the mass fraction in the intimate mixture and particle mean diameter were left as free parameters of the optimisation process, in addition to the macroscopic porosity parameter . For each set of free parameters, a Hapke model was applied to compute the reflectance of the modelled surface, allowing for a direct comparison with New Horizons spectra.
As shown in Fig. 16, the best fits between the modelled surface reflectance and New Horizons spectra were obtained with 1%CH 4 tholins for both H 2 O-rich and H 2 O-poor terrains of Cthulhu (about 44% and 49% of 2 improvement with respect to 5%CH 4 tholins for H 2 O-rich and H 2 O-poor regions, respectively). This is in agreement with the direct comparison between New Horizons spectra and laboratory reflectance measurements of pure tholins discussed in Section 5.2. As already highlighted, the tholins composition (N 2 ∶CH 4 ratio) has a significant influence on the curvature and steepness of the red slope observed in the visible. Modifying the initial composition of the gas mixture from which the tholins are synthesised might therefore be key when trying to reproduce New Horizons spectra. The relative amounts of the refractory species which led to the best fits are reported in Table 2, but should be treated cautiously given the poor quality of the fits and the non-uniqueness of the optimisation solution.
Focusing on H 2 O-rich regions first, the simulated reflectance of a mixing between tholins, hydrocarbon and water ices could not accurately reproduce the New Horizons spectrum, even with 1%CH 4 tholins. More precisely, the MVIC point at = 0.86 μm is not matched and the average reflectance level in the near-infrared of the simulated spectrum is lower than for LEISA spectrum. In particular, this yields a poor fit over the following spectral intervals: 1.3-1.5 μm, 1.65-1.9 μm, and 2.1-2.3 μm, as shown in Fig. 16a. The H 2 O band at 2.0-2.1 μm is also poorly reproduced.
The weighted MVIC data in the visible provide tight constraints on the required tholins ratio. It thus indirectly drives the average reflectance level in the near-infrared. Although the strong H 2 O bands in the near-infrared helps the fitting process by constraining the relative amount of water ice required, the optimisation algorithm fails to reconcile the MVIC data points with high reflectance in the near-infrared.
The fit between the Hapke model and the New Horizons spectrum is significantly worse for the core, H 2 O-depleted regions of Cthulhu (see Fig. 16b). Strong tholins' overtone/combination bands are present in the near-infrared, while they did not appear in the H 2 O-rich spectrum due to the presence of clear H 2 O bands. To reproduce the averaged spectrum collected over H 2 O-poor regions of Cthulhu, H 2 O-ice is indeed no longer needed, but matching MVIC data points requires a significant amount of tholins.
This induces a drop in reflectance around = 1.9 μm, consistent with the pure tholins spectra (see Fig. 12). Around 1.9 μm, the reflectance of the modelled surface thus becomes much lower than for the New Horizons spectrum. In that spectral range, adding more ices to the model would only decrease the reflectance level further because of the hydrocarbon ices absorption bands. Consequently, the optimisation process failed not only to remove the tholins bands which are clearly absent in New Horizons data, but also to reproduce the 2.1 μm and 2.3 μm New Horizons absorption bands (see Section 6.1.3).
In the visible, the difference in the spectral slope is larger at = 0.49 μm, but slightly slower at = 0.86 μm, compared to the H 2 O-rich region case.
It must be noted that scaling the MVIC data down or up by 7% (in agreement with the estimated calibration uncertainty, see Section 2.2) did not yield any noticeable change in the fit between the modelled spectra and New Horizons data.

Ices condensed on tholins particles
With the second surface model described in Section 4.2.3, the fraction of core tholins embedded in the surrounding ice matrix for each icy particle was added to the list of free parameters (already including the mass ratios of each component of the intimate mixing, as well as their respective mean diameters). The best fits obtained by the optimiser with this second surface model are provided in Fig. 17. Fig. 17 compares the fits obtained with the two different surface models (intimate mixing in red, ices condensed around tholins in light blue), for H 2 O-rich and H 2 O-poor Cthulhu regions and with both 5%CH 4 and 1%CH 4 tholins. While no significant difference is noticeable between the two surface models for the H 2 O-poor region case, modelling ices as condensed around core tholins yields a slightly better fit for H 2 O-rich regions (about 10% and 34% of 2 reduction for 5%CH 4 and 1%CH 4 tholins, respectively). However, this limited improvement is to be put in perspective as the discrepancies between New Horizons data and the modelled reflectance spectra remain significant regardless of which model is used, especially for H 2 O-poor regions. It is thus difficult to conclude about how realistic one surface model is compared to another.
Similar conclusions can overall be drawn from the optimisation results, compared to those obtained with the intimate mixing surface model (Section 6.2.1). 1%CH 4 tholins again performed better in reproducing New Horizons data than 5%CH 4 tholins (see Fig. 17). They overall led to an improved fit, especially being closer to MVIC data in the visible and matching the photometric level better in the nearinfrared. The relative mass ratios of the different refractory species corresponding to the best fits for this second surface model are also provided in Table 2.
The reflectance of the H 2 O-rich region of Cthulhu is again reasonably well reproduced (see Figs. 17a and 17b). However, the simulated spectrum still displays strong tholins absorption bands when little water ice is required while they do not appear in LEISA data (i.e. over H 2 Opoor Cthulhu regions, see Figs. 17c and 17d). Moreover, the misfit between the modelled spectra in the visible and MVIC data is far from being fully resolved by using this different surface model, for the point at = 0.86 μm in particular.
It is interesting to note that when using the Maxwell-Garnett model to compute the optical constants of ices condensed around tholins particles, the method selected to derive the value of does not have any significant effect. As detailed in Section 4.2.3, the method chosen in this work is to use the values provided in Khare et al. (1984) and compute the values of which match the single scattering albedo of our laboratory-synthesised tholins. However, the optimisation process was also conducted with a second set of optical constants for the tholins, obtained by assuming remains constant over the spectral range of interest. No noticeable difference was observed between the results achieved with the two separate sets of tholins optical constants.

Discussion
Because of the remaining discrepancies between the modelled surface spectra and New Horizons data (see Section 6.2), we cannot yet draw any firm conclusion about the surface's composition and properties in Cthulhu region, and in particular about the exact chemical nature of the red material. Reliable information about the surface can indeed only be derived from Hapke model's parameters when both the photometric level and spectral features are properly fitted.
Especially, as previously highlighted, the simulated spectra for H 2 Opoor terrains still display tholins absorption features which are absent in LEISA observations. Assuming that the presence of these features is not only due to experimental limitations discussed in Section 6.1.2, three possible explanations are discussed here.
The terrains may be contaminated by interplanetary dust, bringing dark material into the aerosol layer. However, to completely attenuate the tholins bands and thus explain their absence in New Horizons spectra, a significant fraction of darkening material would be needed, which would lead to a severe drop in reflectance (see results presented in Section 5.3). This seems rather inconsistent with the high reflectance observed in the near-infrared in LEISA data. We still tried to reproduce New Horizons spectra using darkened tholins mixed with water and hydrocarbon ices. However, the fit between the modelled surface and the spectra collected by the MVIC and LEISA instruments was worse than the one obtained with pure tholins samples (figures are not shown as the extremely poor quality of the fit does not convey any valuable information). The absence of absorption features thus cannot be explained by interplanetary dust contamination of the red material, as the drop in reflectance it would induce in the near-infrared was shown to be incompatible with the high photometric level measured by LEISA over the Cthulhu region.
As a second hypothesis, Galactic Cosmic Rays (GCR) irradiation is known to promote dehydrogenation reactions, formation of unsaturated bonds (olefinic or acetylenic), cross-linking and disorder (Balanzat et al., 1995). This phenomenon could explain the darkening of tholins and the attenuation of their absorption features, and might therefore justify the absence of absorption bands in New Horizons spectra. This would likely affect the photometric level in the near-infrared as well, and not only the spectral contrast. Radiolysis could also have an effect on the shape of the spectral slope in the visible through the production of conjugated unsaturated carbons, which might alter the current fit of the spectral slope for the H 2 O-poor terrains of Cthulhu. The irradiation dose for a tholins grain on the surface of Pluto can be estimated using the proton flux of Webber and Yushak (1983), assuming the contribution of heavier ions by scaling the proton flux density with their cosmic abundance (Meyer et al., 1998). Here, it is assumed that the GCR density flux is that of local Interstellar Medium, neglecting heliospheric shielding. We also neglect the energy loss of the incoming ions in Pluto's atmosphere. This estimate is therefore an upper value. The electronic and nuclear stopping powers, and , respectively, were calculated with the SRIM software for a H:C:N=1:1:1 stoichiometry and a tholins density = 1 g cm −3 (Ziegler and Manoyan, Table 2 Relative amounts of materials which led to the best fits between the reflectance model and New Horizons spectra, when using 1%CH 4 tholins as an analogue for Cthulhu's dark-red material. The values reported in this table are the mass ratios of all refractory species within the first surface unit (i.e. refractory unit, see Sections 4.2.2 and 4.2.3) and are given in percentage. Results are provided for the two different surface models and for both H 2 O-rich and poor regions. It must be noted that these values must be considered extremely carefully, because of both the non-uniqueness of the solution and poor quality of the fits.
where is the molar mass, the Avogadro number, and the edges of the energy range, the proton flux density from Webber and Yushak (1983), ( ) the relative abundance of element Z relative to hydrogen and is the irradiation time. The irradiation time of a given grain in the aerosol layer depends on the deposition rate and the penetration depth of ions, which is controlled by its mass and energy. Fig. 18 reports the GCR density flux plotted versus the penetration depth for H and Fe ions, which are the lightest and heaviest abundant ions that contribute to the deposited dose. In the case of protons, 99% of GCR energy is deposed within 0-24 m. The aerosol layer (roughly estimated to be of the order of 10 m, Grundy et al., 2018) was then totally irradiated during 4.55 Gyrs, but a given grain within this layer received a varying dose due to the continuous aerosol deposition. In contrast, the iron ions deliver 99% of GCR energy in the range 0-0.48 m, which corresponds roughly to an aerosol deposition duration of 140 Myrs.
The probing depth of visible and near-infrared photons can be estimated to 1 mm, as the mean optical path length calculated from Hapke theory (Clark and Roush, 1984). This depth corresponds to an irradiation time of 280 000 years with the sedimentation rate of Grundy et al. (2018). This leads to upper values of the electronic and nuclear doses deposited in this thin layer of 2.4 × 10 −2 eV atom −1 and 10 −5 eV atom −1 , respectively. The electronic dose is small and may lead to weak or moderate chemical and structural changes, according to previous studies on polymers degradation (e.g. Balanzat et al., 1995;Faure et al., 2021). Therefore, we expect a limited effect on the chemical and optical properties of aerosols. However, these results need to be confirmed by dedicated irradiation experiments on tholins.
Finally, a high porosity of the dark terrains is another possible explanation for the lack of strong absorption features caused by the dark material in New Horizons data. Laboratory experiments showed that a highly porous tholins crust, formed from sublimation experiments, does not display combination/overtone bands in the near-infrared (Poch et al., 2016). It must be stressed that such micro-porosity effects are by definition not accounted for by the macroscopic porosity parameter in Hapke model, and are thus not included in our work. Fig. 19 reports the reflectance spectra of a mixture of water-ice, olivine, N 2 ∶CH 4 =95%:5% tholins (without CO) and smectite, collected before and after sublimation and presented in Poch et al. (2016). The reflectance of our Pluto tholins is displayed as well, and the absence of any remarkable spectral feature in the post-sublimation measurements is clear (from neither tholins nor any of the other species).
Evidences for N 2 and CH 4 sublimation cycles on Pluto's surface were described in Schmitt et al. (2017) and Bertrand et al. (2018Bertrand et al. ( , 2019, but these sublimation cycles do not appear to involve Cthulhu region. This region is actually supposed to be depleted of N 2 -ice (N 2 being too volatile to subsist in these hot, low albedo areas). Moreover, the very small amounts of CH 4 detected in Cthulhu are mainly located on mountains (Bertrand et al., 2020), while there is no spectral evidence for CH 4 in the reddest, H 2 O-poor Cthulhu terrains. Only a very thin CH 4 condensation layer may occur during night (Bertrand et al., 2020).
However, the lack of strong evidence for CH4 sublimation cycles in Cthulhu region after New Horizons' flyby does not completely rule out the possibility of this region being covered by CH4 frosts during other seasons or obliquity periods. Even without sublimation cycles, deposition of light aerosols at very low temperature and under the weak gravity acting on Pluto could anyway promote the formation of highly porous surfaces, which do not show bands in the near-infrared. Again, this hypothesis requires further investigation. Simulating porous surfaces via sublimation of a Pluto tholins/ice mixture or via gentle sedimentation of Pluto tholins at low temperature would be an interesting future step, in order to assess how the porosity affects their reflectance spectrum.

Conclusions
Different hypotheses still coexist to explain the origin and composition of the non-icy, reddish material covering Pluto's dark terrains. They include irradiation of the surface ices (Cruikshank et al., 2015), as well as formation of complex, macromolecular organic compounds in a warm liquid ocean created by the giant impact which formed both Pluto and Charon (Sekine et al., 2017). It has also been suggested that Pluto's dark material could consist of deposited aerosols originating from the dissociation and ionisation of Pluto's atmospheric gases (e.g. Cheng et al., 2017;Grundy et al., 2018). This would be in agreement with the recent findings exposed in Protopapa et al. (2020) according to which a single colouring agent could explain the whole diversity of Pluto's surface colours. In this perspective, this study led to a better characterisation of the reflectance properties of two types of Pluto tholins, to assess their relevance as plausible aerosols analogues for the yet unidentified reddish material covering Cthulhu region.
The reflectance spectra of these tholins were collected under various illumination and observation geometries, highlighting a strong backscattering behaviour. We then numerically retrieved Pluto tholins' single scattering albedo and phase function from their reflectance measurements, by inverting a simplified Hapke model. Feeding these optical properties into a simplified Hapke reflectance model, an optimisation process was set up to assess whether it is possible to reproduce New Horizons spectra with our reflectance model while using Pluto tholins as analogues for the dark material. A mix of tholins and CH 4 , H 2 O, CH 3 OH, and C 2 H 6 ices was assumed and two different surface models were tested (intimate mixing between all refractory species and condensation of icy species around core tholins particles, respectively).
1%CH 4 tholins overall led to significantly better fits than 5%CH 4 tholins. Especially, they better match the MVIC data in the visible and the average photometric level in the near-infrared. This indicates that slightly changing the composition of our Pluto tholins might improve the fit with respect to New Horizons data. Still, our optimisation algorithm failed to fit the three MVIC data points perfectly, even with 1%CH 4 tholins. This misfit in the visible needs to be investigated further, especially since reproducing the strong red slope that Cthulhu's spectra exhibit is key to unravel the composition and origin of the dark red material covering this region of Pluto.
New Horizons data collected over the H 2 O-rich region of Cthulhu are easier to reproduce with our surface reflectance model than H 2 Opoor region' spectra. The strong H 2 O absorption bands dominate the near-infrared part of the H 2 O-rich region' spectrum and hide the absorption bands clearly visible in the tholins spectra. On the other hand, the four absorption features of Pluto tholins appear in the simulated spectra of H 2 O-poor region, while they are absent in New Horizons data.
The presence of the tholins bands in the modelled spectra prevented the optimisation algorithm from reproducing the near-infrared features of the New Horizons spectrum. It is therefore crucial to understand which phenomena could explain such an attenuation of the spectral contrast before concluding about how representative tholins are of the dark red organic material on Pluto. The possible supply in dark material brought by impacting interplanetary dust particles does not seem to be a very plausible hypothesis, since mixing tholins with a darkening agent tends to strongly decrease the reflectance level and not only remove the absorption bands. The effect of GCR irradiation remains to be further investigated, but its contribution might likely be insufficient. However, the effect of a high porosity in the Cthulhu terrains might provide an explanation for the lack of strong absorption features in New Horizons spectra, and also requires deeper analyses.
The presence of some of the absorption bands in Pluto tholins' reflectance spectra could also be partially related to the way they are synthesised. This more precisely concerns the two bands related to amine functions (i.e. 1.5 μm and 2.0 μm bands, see Section 6.1.2). The formation of amine functional groups is indeed known to be enhanced by wall effect in the PAMPRE experimental setup, which could explain the intensity of some of the absorption bands. Another point is that the tholins synthesis was conducted at room temperature, while the formation of haze particles occurs at much lower temperature in Pluto's atmosphere. To our knowledge, no tholins were collected at low-temperature in cold plasma reactors. Bernard et al. (2003) and Bernard et al. (2006) and Brucato et al. (2010) used a cryogenic trap in their experiments to collect volatile species, but their samples were synthesised at room temperature and cannot be used as low-temperature tholins (as mentioned in Brasse et al., 2015). Some temperature effects are however expected, such as the condensation of free molecules that would result in an increase of the soluble organic fraction. Nevertheless, we still expect the presence of similar functional groups as in tholins synthesised so far and we therefore do not foresee any significant impact on our conclusions.
Finally, all our analysis is based on extremely limited reflectance models, the simplicity of which might affect our simulated spectra. The isotropic phase function assumption may for instance alter the photometric level, but also the spectral slope and the band depths. Pluto tholins show a strong backscattering behaviour (see Section 5.2), while ices are on the contrary characterised by a strong forward-scattering. The variations of the phase function as a function of both the geometry and the wavelength are thus hard to predict, making the isotropic assumption a severe limitation of our current model.
It must also be stressed that there is no guarantee of the solution's uniqueness when fitting a Hapke reflectance model to an existing spectrum. Nonetheless, this does not hinder the relevance of our approach since we are not estimating specific surface properties from our best fit, but are only trying to reproduce Cthulhu's reflectance using tholins as analogues for Pluto's dark organic material.
Finally, the optical properties of our Pluto tholins (including optical constants and ) were retrieved from a purely numerical inversion of a highly simplified radiative transfer model, applied on reflectance measurements performed on a complex medium. However, to rigorously extract optical constants, it is preferable to use films deposited on plane surfaces, for which simplified optical transmission models are bettersuited (Trotta, 1996). For this reason, the inverted optical constants used in this work and provided in Appendix A.3 need to be treated cautiously. An additional, similar analysis based on recently determined optical constants for Pluto tholins might be insightful (Jovanović et al., 2021).
This appendix describes how to invert Hapke model, in order to extract some surface parameters of interest from reflectance measurements.

A.1. Estimating the single scattering albedo and phase function
As mentioned in Section 5.5, a linearised least-squares method can be used to invert Hapke model and retrieve the single scattering albedo and phase function from the tholins reflectance spectra.
The linearised least-squares method is by definition an iterative process. The non-linear relation between the parameters to estimate (parameters vector ) and the observations (here reflectance values, concatenated in an observations vector ) can be written as follows: where the matrix ( ) represents the available model linking the observations to the parameters. From the Taylor series of the above equation, a first order approximation can be derived for the deviation in the observations which is expected from a small change in the parameters: where ⋆ is the current best estimate of the parameters vector . To simplify the notations, the partial derivatives matrix (also referred to as the normal matrix) is introduced, as: The previous step (Eq. (15)) is necessary to linearise the problem. The successive updates which need to be applied to the parameters vector, starting from an initial guess 0 , are given by where is the difference between the actual reflectance data and the simulated reflectance values obtained with the former estimate of the parameters vector .
The parameters to be estimated in our particular problem are the single scattering albedo (as a function of the wavelength ) and the phase function which depends on the phase angle . The model relating ( ) and ( ) to the expected reflectance is the simplified Hapke model (see Eq. (1)). The least-squares inversion as described in Eq. (17) requires the computation of the normal matrix defined by Eq. (16). Consequently, the partial derivatives of the simulated reflectance ( , , , ) with respect to both ( ) and ( ) need to be computed. Let us first define the function ( , 0 , ) as the product of the two multiple scattering terms in Eq. (1): Then, the derivative of the reflectance with respect to the single scattering albedo, for a given value of the wavelength , can be computed as follows: where the partial derivative of the function ( , 0 , ) can be expressed as Finally, the partial derivative of the reflectance with respect to the phase function is assuming the phase angle is constant. After a few iterations, the linearised least-squares inversion described above converges towards an estimated solution for the single scattering albedo and the phase function.

A.2. Extracting the optical constants
As already discussed in Section 4.2.3, the Maxwell-Garnett model requires the optical constants of the tholins as input. However, the least-squares inversion of the simplified Hapke model described in Appendix A.1 only provides an estimate of the single scattering albedo ( ) and phase function ( ) for the laboratory synthesised tholins. The highly non-linear relation between the single scattering albedo and optical constants and (see Eq. (2)) does not allow for a straight-forward determination of the tholins optical constants from .
The spectral information is mostly carried by , except for strong fundamental bands. The near-infrared spectral range covered by LEISA values for both 5%CH 4 and 1%CH 4 tholins, obtained from the numerical inversion described in Appendix A.2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) spectra corresponds to the domain of overtone and combination bands. Variations of the refractive index are therefore expected to not exceed a few percents. On the other hand, the very strong band in the UVvisible responsible for the red slope should lead to a large change in the refractive index in the visible. Using the values provided in Khare et al. (1984) aims at taking this variation into account.
From these assumed values, we estimate the corresponding values of that would match the single scattering albedo . To compute as a function of the wavelength, the expression giving the single scattering albedo as a function of and (Eqs. (7) and (8)) must be inverted for each wavelength. Because the expression for is highly non-linear, a linearised least-squares method can be employed, exactly similar to the one presented in Appendix A.1.
The observations are here the single scattering albedo values for all wavelengths of interest, while the parameters to be estimated are the associated values. The normal matrix contains the partial derivatives of with respect to , evaluated separately at each wavelength. The iterative linearised least-squares method converges within a few iterations and provides estimates for the values.

A.3. Tabulated results of the Hapke model inversion
This appendix provides the tabulated values for the single scattering albedo and phase function obtained after performing the inversion of the Hapke model as described in Appendix A.1. Results are reported for both 5%CH 4 and 1%CH 4 tholins, in Tables A.1 and A.2 for and , respectively. Table A.3 contains the optical constants extracted from the single scattering albedo, as described in Appendix A.2. values are also displayed as a function of the wavelength in Fig. A.1.
Again, it must be highlighted that those parameters are determined in a purely numerical way. They need to be used carefully, keeping in mind the inherent lack of physical meaning behind the least-squares inversion of our models.

Appendix B. Optimisation process
This appendix provides more details about the optimisation problem introduced in Section 4.3 and about the optimisation algorithm selected to tackle it.

B.1. Definition of the optimisation problem
For this particular problem, the optimisation objective (i.e. variable to be minimised) was defined as the 2 between the simulated reflectance and the New Horizons spectrum.
As a reminder, below are listed the design parameters (i.e. free parameters) whose values were tuned by the optimisation algorithm to find the best match between the simulated spectra and the New Horizons data: • mass ratio of each of the chemical species within the intimate mixture of the refractories spatial unit (H 2 O, C 2 H 6 , CH 3 OH, tholins); • mean particle diameter for each of the chemical species (H 2 O, C 2 H 6 , CH 3 OH, tholins, CH 4 ); • spatial ratio between the two units (the refractories and the volatiles ones); • spatial fraction of tholins in icy grains when the Maxwell-Garnett model is used (Section 6.2.2). If so, the mean particle diameter and mass fraction associated with tholins are removed from the list of the design variables; • macroscopic porosity parameter .
All the mass and spatial ratios are by definition ranging from 0 to 1, which is also the case for the porosity parameter . Additionally, the sum of the mass fractions of the refractory species covering one of two spatial units (see Figs. 2 and 3) must be equal to 1. Finally, the mean particle diameters were constrained between 0.1 μm and 100 μm.

B.2. Selection of the optimisation algorithm
To conduct the optimisation process, we used the PAGMO toolbox developed by the Advanced Concepts Team (ACT) at ESA, which offers a wide range of optimisation algorithms (in C++). The documentation can be accessed at https://esa.github.io/pagmo2/.
As mentioned in Section 4.3, the performance of several global optimisation algorithms was assessed in terms of convergence and computational efficiency for our particular optimisation problem. The simple genetic algorithm, differential algorithm, predator-swarm algorithm, and artificial bee colony were tested. Those are all population-based evolutionary algorithms.
Except for the artificial bee colony algorithm, the other ones all achieved convergence (verified by reconducting the optimisation process with different random seeds). The simple genetic algorithm offered the best trade-off: it seemed to provide a good approximation of the global optimum within a reasonable computational time. The efficiency of the genetic algorithm was of particular relevance for this work, as both the number of design parameters and the associated search space were large, while computational capabilities were limited. For these reasons, this algorithm was selected to tackle the optimisation problem described in Appendix B.1.
The simple genetic algorithm relies on Darwin's natural selection principle, and more precisely on the ''survival of the fittest'' concept (e.g. Goldberg, 1989). As for any other evolutionary algorithm, an initial population of individuals (i.e. solutions, or combinations of design parameters) is generated. The ''fitness'' of each of these individuals is computed, as being equal to the value of the optimisation objective corresponding to the particular set of design parameters defining the individual. The population then evolves to form a second generation. Mimicking the evolution of populations according to evolutionary theories, parents are selected among the current population of solutions and each child individual is computed as a mix between the two sets of design parameters originating from its two parents. Following the natural selection principle, individuals associated with a good fitness value (i.e. for which the value of the optimisation objective is lower) are more likely to be selected as parents among a given generation, and thus to propagate the good performing design parameters defining them to the next generation. As generations keep being generated, the best identified fitness value tends to get closer to the global optimum and the ''fittest'' individuals thus approach the best set of design parameters. To guarantee a good sampling of the research space, the population size (i.e. number of individuals per generation) should not be chosen too small. Additionally, mutations are introduced from one generation to the next to increase the diversity within the population, again based on what is observed in natural species evolution.