# Impacts of Organic Ice Condensation on the Optical Properties of Haze on Pluto

Jingyu Wang, Siteng Fan, Chao Liu, Vijay Natraj, Leslie A. Young, Yuk L. Yung

► To cite this version:

Jingyu Wang, Siteng Fan, Chao Liu, Vijay Natraj, Leslie A. Young, et al. Impacts of Organic Ice Condensation on the Optical Properties of Haze on Pluto. The Planetary Science Journal, 2023, 4, 10.3847/PSJ/acaf30. insu-03993949

HAL Id: insu-03993949

https://insu.hal.science/insu-03993949

Submitted on 17 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
OPEN ACCESS

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

CrossMark

# Impacts of Organic Ice Condensation on the Optical Properties of Haze on Pluto

Jingyu Wang<fnref n="1" />,<fnref n="2" /> ID, Siteng Fan<fnref n="3" />,<fnref n="4" /> ID, Chao Liu<fnref n="1" /> ID, Vijay Natraj<fnref n="5" /> ID, Leslie A. Young<fnref n="6" /> ID, and Yuk L. Yung<fnref n="4" />,<fnref n="5" /> ID

Received 2022 September 29; revised 2022 December 14; accepted 2022 December 22; published 2023 January 27

## Abstract

The flyby of the New Horizons spacecraft in 2015 July revealed an unexpected cold atmosphere of Pluto and confirmed the existence of its atmospheric haze. The observed and simulated vertical profiles of chemical species and microphysical processes suggest that the haze particles in Pluto's middle and lower atmosphere may contain organic ice condensation. Such organic ice components can potentially affect Pluto's haze chemistry and optical properties, as well as its energy budget. This study investigates the influence of the ice components on the scattering properties of Pluto's haze by comparing New Horizons observations and simulated particle scattering properties. Comprehensive tests are performed for various haze particle parameters, including their size, chemical component, ice content, and morphology. Scattering properties of these ice-bearing haze particles are calculated by a discrete dipole approximation method and compared to multispectral observations obtained by four New Horizons instruments in spectral regions ranging from the ultraviolet to the near-infrared. The results indicate that the inclusion of the organic ice component leads to higher ratios of backscattering in the visible to extinction in the ultraviolet and provides better agreement with observations compared to monodispersed homogeneous aggregates. But it alone is not sufficient to explain the observed forward scattering values in the visible and near-infrared. Therefore, other scattering sources and/or mechanisms are still required to explain the full set of scattering observations. Further observations, as well as laboratory measurements and numerical tests, are anticipated to improve our understanding of the morphology and ice content of Pluto's haze.

*Unified Astronomy Thesaurus concepts: Atmospheric composition (2120); Planetary atmospheres (1244); Pluto (1267)*

## 1. Introduction

In 2015, the New Horizons spacecraft's flyby of Pluto measured its scattered radiation from the ultraviolet (UV) to the near-infrared (NIR) and confirmed the existence of haze in the atmosphere (Stern et al. 2015; Gladstone et al. 2016). Fundamental microphysical properties of these particles are then constrained using extensive information of extinction and scattering intensities from multiwavelength and full-phase observations obtained by four instruments on board the New Horizons spacecraft: the Multispectral Visible Imaging Camera (MVIC; Kutsop et al. 2021; Fan et al. 2022), the LOng-Range Reconnaissance Imager (LORRI; Cheng et al. 2017), the Linear Etalon Imaging Spectral Array (LEISA; Grundy et al. 2018), and the Alice UV Imaging Spectrograph (Young et al. 2018).

Pluto's haze extends from the surface to several hundred kilometers in altitude. The bluish color obtained by MVIC indicates the dominance of Rayleigh scattering by particles much smaller than visible wavelengths, while the large high-to-low-phase brightness ratio suggests stronger forward scattering, and thus the existence of large particles (Gladstone et al. 2016). Given these two pieces of evidence, Pluto's haze particles are suggested to be composed of aggregates, i.e., fluffy, and porous

particles consisting of small subunits or so-called "monomers." Homogeneous haze particles that have the same size and composition are often assumed in order to derive their morphology. Gao et al. (2017) explained the observed UV extinction profiles using fractal aggregates resulting from microphysical procedures such as coagulation and transport. However, discrepancies arise when full-phase scattering intensities from LORRI are considered (Cheng et al. 2017). Kutsop et al. (2021) combined the UV extinction obtained by Alice and the scattering intensities at three wavelengths and eight phase angles obtained by MVIC and proposed that the haze size distribution follows either the bimodal or the power-law distribution. By testing a number of scenarios of homogeneous particles' size distributions of Pluto's haze, Fan et al. (2022) confirmed that a bimodal distribution of "tholin"-like particles, including a large population of smaller-sized spheres with a radius of ~80 nm and a small population of larger-sized fractal aggregates with an effective radius of ~1 µm, could result in scattering properties that agree with the observations obtained by the four New Horizons instruments.

With Pluto's atmospheric temperature as low as ~70 K (Wong et al. 2017; Gao et al. 2017; Lavvas et al. 2021), heterogeneous condensation of organic gases is expected to occur in Pluto's atmosphere, with the "tholin"-like aggregates serving as condensation cores. Haze particles generated at higher altitude and transported downward can serve as the condensation cores. As a result of photochemical processes

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters/China Meteorological Administration Aerosol-Cloud-Precipitation Key Laboratory, School of Atmospheric Physics, Nanjing University of Information Science & Technology, Nanjing 210044, People's Republic of China
chao_liu@nuist.edu.cn</footnote>

<footnote n="2">² ETH Zurich, Zurich 8094, Switzerland</footnote>

<footnote n="3">³ LMD/IPSL, Sorbonne Université, PSL Research Université, École Normale Supérieure, École Polytechnique, CNRS, Paris F-75005, France</footnote>

<footnote n="4">⁴ Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="5">⁵ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA</footnote>

<footnote n="6">⁶ Southwest Research Inst., Boulder, CO 80302, USA</footnote>
involving N₂ and CH₄, it is possible to generate “tholin”-like aggregates at altitudes above 400 km (Khare et al. 1984; Gladstone et al. 2016; Wong et al. 2017; Young et al. 2018). The observed local density profiles of organic gases (C₂H₂, C₂H₄, and C₂H₆) reach local minima near 200 km altitude (Young et al. 2018), which cannot be explained by the homogeneous chemistry simulations and instead is attributed to heterogeneous processes on the haze particles. Therefore, organic gases, such as C₂ species and C₄H₂, may condense on them in the middle atmosphere (~200 km), and considerable organic ice components are expected to be mixed with “tholin”-like aggregates (Wong et al. 2017; Young et al. 2018). Lavvas et al. (2021) used condensed C₄H₂ aggregates as the proxy, computed corresponding optical properties using the particle scattering model developed by Rannou et al. (1999), and found better agreement with observations than if homogeneous “tholin”-like particles were used. However, considerable differences still exist in the backward scattering near the surface. Fan et al. (2022) presented a sensitivity test by substituting the “tholin”-like component with organic ice and proposed that such replacement does not significantly influence their scattering properties. However, the aforementioned efforts either qualitatively discussed or oversimplified the influence of gas condensation on the haze scattering properties, leaving the characteristics of organic ice on Pluto’s haze particles still an open question.

Part of the lack of understanding may be caused by inaccurate simulation of optical properties. In simulating scattering properties of aggregates, Rannou et al. (1999) developed a semiempirical model based on Mie scattering and empirical rules, and Tomasko et al. (2008) used a parameterization method based on a set of T-matrix calculations, and these two classic models will be referred to as “R99” and “T08” for simplicity in this work. These two models contain oversimplified assumptions and are not capable of simulating ice-bearing particles. Meanwhile, a number of numerical models have also been developed in the Earth and planetary science community to simulate scattering properties of nonspherical or inhomogeneous particles, such as the T-matrix method (Waterman 1965), Finite-Difference Time-Domain (FDTD) method (Yee 1966), Multiple Sphere T-Matrix method (MSTM) model (Mackowski & Mishchenko 1996), Pseudo-Spectral Time-Domain (PSTD) method (Liu 1997), and the invariant imbedding T-matrix method (Bi et al. 2013). These state-of-the-art numerical models make more accurate simulations of sophisticated Pluto’s haze particles possible. This study uses a discrete dipole approximation (DDA) method that has the capability to simulate the optical properties of particles of Pluto’s haze with inhomogeneous structures. The DDA is a numerically rigorous method to compute particles’ scattering properties by solving Maxwell’s equations (Purcell & Pennypacker 1973; Yung 1978; Draine 2003; Yurkin & Hoekstra 2007; Teng et al. 2019), where particles are treated as dipoles much smaller than the incident wavelength. Thus, it can consider particles with complex morphology and mixing states and is a powerful tool to investigate the influence of gas condensation on Pluto’s haze particle scattering properties.

This study will evaluate influences of the condensed organic ice on the optical properties of Pluto’s haze using the New Horizons observations. The remainder of the paper is arranged as follows. Section 2 introduces the methodologies and

observations. Systematical sensitivity tests for various haze particle and coating component properties are presented in Section 3. The model results are compared with the observations to provide the best-fit scenario in Section 4. Section 5 discusses the implications of the results.

## 2. Methodologies and Observations

This study generates numerical particles by defining their detailed size, morphology, and components. Once a particular particle is built, their corresponding scattering properties will be simulated and compared with the observations. This section details the particles we build, the numerical simulations, and the observations.

### 2.1. Haze Morphology

Following previous studies, “tholin”-like aggregates are supposed to form as a result of photochemistry and coagulation in the upper atmosphere, and then the organic gases condense on them as ice coating when the particles are transported to gas condensation regions in the middle atmosphere near ~200 km (Lavvas et al. 2021). Thus, ice-coated haze particles are expected in the lower atmosphere at <50 km. Considering the vertical profiles of temperature and chemical species abundances, the ice coatings are possibly C₂H₂, C₂H₄, or C₄H₂, or a combination of them (Wong et al. 2017; Lavvas et al. 2021). Therefore, the haze particles in this study will be assumed to be inhomogeneous with a “tholin”-like aggregate coated by organic ice sphere.

The observed bluish color and strong forward scattering indicate that Pluto’s haze mainly consists of fractal aggregates formed through coagulation of small spherical monomers (Gladstone et al. 2016). There are three key variables to parameterize the morphology of such aggregates: fractal dimension ($D_f$), number of monomers in each aggregate ($N_m$), and the radius of the monomer ($R_m$). $D_f$ describes the aggregate porosity. Particles with smaller $D_f$ are fluffier and more porous. This quantity is defined through the mass-radius relation (Gmachowski 2002). The fractal dimension is the exponent of the power and can be determined from the slope of the log-log plot.

Then, an aggregate volume effective radius ($R_a$) is given by

$$R_a = \sqrt[3]{N_m} R_m. \qquad (1)$$

The aggregate materials are physically represented by their complex refractive indices, i.e., $RI_{agg} = N_{agg} + i N_{agg}$.

For simplification, the organic coating is considered as sphere enclosing the aggregate and has the same center as that of the aggregate. In this way, the inhomogeneous particle is simply represented, and two more parameters are sufficient to specify the particle. Coating amount will be given by the volume fraction of the aggregate ($V_p$), and the refractive index of the coating material is represented by $RI_{coat}$. The organic ice coating is assumed to be nonabsorbing, as organic ices have weak absorptivity in the visible and NIR wavelengths, with imaginary index <10⁻⁶ (Khare et al. 1989) and have a representative real refractive index corresponding to C₄H₂ ice, which is the most abundant organic ice in Pluto’s atmosphere as suggested by Lavvas et al. (2021), i.e., $RI_{coat} = 1.4$ (Khanna et al. 1988; Haynes 2013).

Aggregates are numerically generated using the model described in Liu et al. (2019), and concentric spheres are
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Four examples of numerically generated coated haze particles considered in this study. (a) $R_m = 15$ nm, $N_m = 300$, $D_f = 2.0$, $V_f = 0.5$; (b) $R_m = 20$ nm, $N_m = 600$, $D_f = 1.8$, $V_f = 0.9$; (c) $R_m = 15$ nm, $N_m = 880$, $D_f = 2.0$, $V_f = 0.8$; (d) $R_m = 10$ nm, $N_m = 400$, $D_f = 2.5$, $V_f = 0.7$. The small black spheres denote monomers, and the blue sphere denotes the ice coating.</figcaption></figure>

added. Four examples of numerically generated coated aggregates with different aggregate morphologies and coating properties are shown in Figure 1. Notice that the ways organic ice condenses on the aggregate can be complex. There may be condensation on each monomer (Teng et al. 2019), or the coating can be irregular (Liu et al. 2015; Dong et al. 2015). In this work, we use the simple coating model to investigate the influence of the presence of coating.

Most haze particle properties mentioned above are tested. As a benchmark, we adapt parameters of uncoated aggregates from Fan et al. (2022), the best constraints on Pluto's haze available to date. As forward scattering is sensitive to the particle size parameter, we pick (for the initial guess) a set of parameters from Fan et al. (2022) retrieved by only considering forward scattering; we then conduct further tests to evaluate the contribution of the ice coating. The parameters $D_f$, $R_m$, and $N_m$, are assumed to be 1.9, 24 nm, and 1275, respectively. The scattering properties of such haze particles are first computed using the DDA method. Then, several combinations of different particle parameters are examined to investigate their influences on the optical properties of the inhomogeneous haze particles.

A summary of the values and ranges of the tested parameters is provided in Table 1. The value of $D_f$ varies from 1.8 to 2.8. $D_f$ of freshly formed aggregates is expected to be around 1.8 in Earth's atmosphere, and aging processes will make the structure more compact, resulting in a typical value of around 2.8 (Li et al. 2016; Wu et al. 2016). The value of $R_m$ ranges from 10 to 30 nm, a typical monomer size for Pluto's haze (Gao et al. 2017). Seven volume fractions of the aggregate from 0.25 to 1.00 are considered, with $V_f = 1$ representing uncoated

<table><caption>**Table 1**
Summary of the Values of the Parameters Tested</caption><thead><tr><th></th><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td rowspan="4">Aggregate</td><td>$D_f$</td><td>1.8, 2.0, 2.2, 2.5, 2.8</td></tr><tr><td>$R_m$ (nm)</td><td>10, 15, 20, 25, 30</td></tr><tr><td>$V_f$</td><td>1, 0.90, 0.80, 0.70, 0.60, 0.50, 0.25</td></tr><tr><td>$RI_{agg}$</td><td>Khare et al. (1984); Jovanović et al. (2021)</td></tr><tr><td rowspan="2">Coating</td><td>$RI_{coat}$</td><td>1.3, 1.4, 1.5, 1.6</td></tr><tr><td>Shape</td><td>sphere</td></tr></tbody></table>

aggregates. Two sets of data for the refractive indices of aggregates are considered. The first is widely used in simulating Titan's haze and is adapted to the Pluto context as well (Khare et al. 1984). The other is from a recent laboratory measurement using a setup representative of Pluto's atmosphere (Jovanović et al. 2021). The real parts of the refractive indices of the two are similar, while Jovanović et al. (2021) suggests larger imaginary parts at wavelengths smaller than 380 nm and larger than 900 nm, likely due to the higher oxygen fraction in Pluto's haze particles. The optical properties from Jovanović et al. (2021) are limited to wavelengths from 270 to 2100 nm, so the refractive index at 185 nm is taken from Khare et al. (1984). For the nonabsorptive coating, real parts of the refractive index ($N_{rcoat}$) ranging from 1.3 to 1.6 are tested, which represent condensation of different organic gases.

## 2.2. The DDA Method

The DDA method serves as a rigorous model to compute the optical properties of complex inhomogeneous particles,
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Ratios of scattering intensities to UV extinction at the wave bands of (a) LORRI, (b) LEISA, and (c–f) MVIC for different aggregate fractal dimensions ($D_f$). Results from simulation of different aggregate fractal dimension values are denoted by colored lines; the black curves indicate results for the reference particle derived from Fan et al. (2022), which is used as the initial guess in this work.</figcaption></figure>

<table><caption>**Table 2**
Observations of Pluto's Haze</caption><thead><tr><th>Instrument</th><th>Wavelength (μm)</th><th>Altitude (km)</th><th>Phase Angle (deg)</th><th>Data</th><th>Number of Observations Considered in This Work</th><th>Reference</th></tr></thead><tbody><tr><td>Alice</td><td>0.185</td><td>0–300</td><td>0 (Extinction)</td><td>Extinction coefficient (α<sub>ext,UV</sub>)</td><td>18</td><td>Young et al. (2018)</td></tr><tr><td>LORRI</td><td>0.608</td><td>0–100</td><td>19.5</td><td>Local scattering intensity (I<sub>sca,vis</sub>)</td><td>36</td><td>Cheng et al. (2017)</td></tr><tr><td></td><td></td><td>0–50</td><td>67.3</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>0–75</td><td>148.3</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>0–200</td><td>169.0</td><td></td><td></td><td></td></tr><tr><td>LEISA</td><td>1.235–2.435</td><td>0–299</td><td>169.0</td><td>Local scattering intensity (I<sub>sca,NIR</sub>)</td><td>225</td><td>Grundy et al. (2018)</td></tr><tr><td>MVIC</td><td>0.492 0.624<br>0.863 0.881</td><td>0–50</td><td>18.2 38.8 169.4</td><td>Local scattering intensity (I<sub>sca,MVIC</sub>)</td><td>108</td><td>Fan et al. (2022)</td></tr></tbody></table>

including the organic-ice-coated haze particles employed in this study. In this model, particles are divided into cubic subvolumes, which are called dipoles. Interactions between dipoles are computed using integral equations for the electric field, and then the polarization of the scatter is solved, which leads to the final scattering quantities. Optical properties given by DDA include scattering phase functions, cross sections, and efficiencies of particles at the relevant wavelengths. The scaling

factor in the model, the number of dipoles per wavelength (dpl), is assumed to be three times the ratio of wavelength to monomer radius; this represents a compromise between accurately characterizing the particle geometry and the computational cost. For instance, a typical range for dpl for a particle with a radius of 15 nm in the NIR (1.235–2.435 μm) is from 247 to 487. Details of this model can be found in Yurkin & Hoekstra (2007).
<figure><img src="image_6.png" /><figcaption>Figure 3. Same as Figure 2, but for different monomer radii R<sub>m</sub>.</figcaption></figure>

Several models have been developed to simulate scattering properties of nonspherical particles. However, compared to the R99 and T08 models that are widely used for planetary haze particles, DDA is capable of simulating inhomogeneous and arbitrarily shaped particles, such as coated aggregates considered in this work. Since particles are discretized into dipoles in the model, it can represent the aggregate and coating separately and thus can treat particle morphology and components accurately without oversimplification. Methodologically, results from DDA are expected to be more accurate because DDA directly solves Maxwell’s equations in the frequency domain, while R99 and T08 oversimplify the scattering process and employ assumptions such as the omission of polarization (Rannou et al. 1999). Furthermore, for aggregates with more than 10³ monomers, DDA is more reliable than R99 and T08, whose maximum test limits are 1000 and 1024, respectively. Therefore, we choose to employ the DDA method in this study to accurately represent the coated aggregate.

## 2.3. Observations

To better constrain the morphology and component of haze particles with ice components, we consider observations from four instruments on board the New Horizons spacecraft (Fan et al. 2022). A summary of the observations, including the observational wavelength, phase angle, and altitude range, is

given in Table 2. Among the four instruments, Alice and LEISA provide spectral observations in the UV and NIR regions, respectively, while LORRI and MVIC are broadband and narrowband cameras, respectively, measuring in the visible. Notice that the Alice measurements are correlated, since the analysis of Young et al. (2018) smoothed transmission in altitude before retrieving haze properties. Because the error in the ratio is dominated by the numerator (local scattering intensity), this has little practical effect on the analysis.

The scattering features of Pluto’s haze in the 5–50 km altitude range are similar (Fan et al. 2022), so all observations within the range are used together to evaluate the parameter combinations. The agreement between observations and modeled results is quantified by a single variable χ², defined as

$$ \chi^2 = \sum_{i=1}^{n} \sum_{j=1}^{m} \left[ (RS_{ij} - R_{ij})^2 / \sigma_{D_{ij}}^2 \right], \quad (2) $$

where i and n are the index and total number of altitude bins, respectively; j and m are for independent observation variables. Here RS<sub>ij</sub> and R<sub>ij</sub> are, respectively, the simulated and observed ratio of scattering intensity to UV extinction, and σ<sub>D<sub>ij</sub></sub> is the observation uncertainty at the ith altitude for the jth observation. The altitude bins are centered from 7.5 to 47.5 km with an interval of 5 km, resulting in nine bins. The number of
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Figure 2, but for the ratio of scattering intensity to UV extinction at varying aggregate refractive index (RI<sub>agg</sub>).</figcaption></figure>

observations from each instrument is shown in Table 2; there
are a total of 387 observations considered for the evaluation.

**3. Influence of Individual Parameters**

Systematic sensitivity tests are conducted to investigate the
influence of particle morphology and composition (R<sub>m</sub>, D<sub>f</sub>, V<sub>f</sub>,
RI<sub>agg</sub>, and RI<sub>coat</sub>) on their scattering properties. For each
parameter varied in the sensitivity test, other parameters remain
the same as our initial guess mentioned in Section 2.1 (D<sub>f</sub>, R<sub>m</sub>,
and N<sub>m</sub>, are 1.9, 24 nm, and 1275, respectively; note also that
the initial guess assumes no coating). The simulated results are
presented as ratios of scattering intensity to UV extinction at
185 nm. The ratio is defined as R = (1/4)P<sub>11</sub>C<sub>sca</sub>/C<sub>ext</sub>. We use
the average values of the UV extinction obtained by Alice
during solar egress and ingress, which differ by less than 10%.
The observational uncertainties based on error propagation are
also considered.

Figure 2 illustrates the influence of aggregate fractal
dimension D_f. With a larger D_f, aggregates become more
compact, and their forward scattering intensities become
stronger at all wave bands, with a slight increase in the
intensities at backward scattering angles. However, the slope of
the phase function curve in the forward scattering gets smaller.
Thus, with the aging of the haze particles as they get
transported, they tend to have larger forward scattering
intensities.

Similarly, the influence of R<sub>m</sub> is shown in Figure 3. Larger
monomers lead to stronger scattering intensities at all scattering
angles for a given UV extinction; this holds for both the visible
and infrared wavelengths. The intensity of forward scattering
becomes larger as compared to backward scattering, resulting
in more asymmetrical phase functions that are more evident in
the LORRI wave bands (Figure 3(a)) than those of MVIC
(Figures 3(c)–(f)). This indicates that when the intensity of the
phase function to be fitted is larger than the assumed value at
all scattering angles, or when the phase function curve is more
asymmetrical, increasing the monomer radius can improve
the fit.

The aggregate and coating components, in the form of their
RI, are also essential factors determining particle optical
properties. Figure 4 shows the scattering properties of
"tholin"-like particles based on the refractive indices given
by Khare et al. (1984) and Jovanović et al. (2021). The RI
measurements have very similar values of N_r around 1.65,
while the N_i values found by Jovanović et al. (2021) are one
order of magnitude higher in the NIR. As the UV measure-
ments are not included in Jovanović et al. (2021), the UV
extinction for both assumed analogs is computed using the data
from Khare et al. (1984). Nevertheless, the differences in the
haze particle scattering properties due to the two sets of RI<sub>agg</sub>
are negligible. This indicates that the differences in RI<sub>agg</sub>
presented by the two data sets are not sufficient to have
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Figure 2, but for the ratio of scattering intensity to UV extinction at varying coating refractive index (RI<sub>coat</sub>). The volume fraction of the aggregate is fixed at 0.75.</figcaption></figure>

noticeable influence on the New Horizons observations, which agrees with previous studies (Fan et al. 2022).

Figure 5 shows the influences of RI<sub>coat</sub>. The chemical composition of the organic components is suggested to be C<sub>2</sub> species or C<sub>4</sub>H<sub>2</sub> (Lavvas et al. 2021). As the exact RI of organic ice is not well reported, a range of RIs (1.3, 1.4, 1.5, and 1.6) around values representative of typical organic ices is considered. Results show that the larger the real part of the RI, the stronger the scattering intensity of the ice coating, which causes the altitude-averaged difference to increase. This is seen in the visible and NIR results. Therefore, aside from increasing the monomer radius, changing RI<sub>coat</sub> of the coated aggregates can also alter the scattering intensity at every scattering angle, but with a smaller influence on the asymmetry factor.

The volume fraction of the aggregate ($V_f$) is the ratio of the volume of the aggregate to that of the haze particle, so particles with smaller $V_f$ contain more ice coating. Two tests are performed and presented in Figures 6 and 7. Figure 6 shows results for particles with a fixed aggregate and with increasing ice coating amount, i.e., $V_f$ decreasing from 1.0 to 0.25. A larger amount of ice coating increases the scattering intensity but shows less influence on the shape of the scattering phase function curve. Results shown in Figure 7 are for coated particles that are formed by replacing part of the aggregate with organic ice while keeping the particle size parameter fixed. When the particle size parameter is fixed, with increasing coating volume fraction, the slope of the curve in the forward

scattering range decreases, while that in the backward scattering range increases slightly. Ice coating leads to stronger scattering at all angles and wavelengths. Thus, changing $V_f$ with fixed particle size parameter will also alter the asymmetry factor, while adding coating directly onto the aggregate can lead to larger scattering intensities, which is similar to the effect of increasing RI<sub>coat</sub>.

Therefore, ice coating can influence the shape of the scattering phase curves. Previous studies assumed various size distributions of haze particles, which needs more evidence based on the haze production mechanism. More importantly, inhomogeneous particles and the morphology with particular coating and core were less discussed, making the impacts of the ice coating on the particle's scattering properties less known. In order to investigate the influence of ice coating, as well as how well the coating may improve the agreement between simulated and observed properties when employing a simpler assumption of one population of coated aggregates, we use grid search to further constrain the parameters describing the particle.

4. Best-fit Solution

To examine the extent to which coated haze particles can explain the observations that potentially include implications on Pluto's haze properties, we conduct a series of grid searches among the five key parameters considered in Section 3 ($N_m$, $R_m$, $D_f$, $V_f$, and RI<sub>coat</sub>). Values of RI<sub>agg</sub> are from Khare et al. (1984),
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Figure 2, but for the ratio of scattering intensity to UV extinction at varying volume fraction of the aggregate ($V_f$) with a spherical coating for simple addition.</figcaption></figure>

with the real part ranging from 1.62 to 1.71 and the imaginary part ranging from $4 \times 10^{-4}$ to 0.22. The searches are performed alternatively with two types, first varying the parameters that describe uncoated aggregates, and second adding coating onto the aggregates, in order to obtain the best fit. Two additional parameters describing the coating ($V_f$ and RI<sub>coat</sub>) are investigated to generate the coated haze particle that fits best with observations.

First, for the simulation of uncoated aggregates, we consider the observations of forward scattering in visible and IR and extinction in UV. The three parameters describing the particle's morphology are $N_m$, $D_f$, and $R_m$. The ranges of these three parameters are smaller than those used for the sensitivity test (Section 2, Table 1) based on the result from it. In the first iteration of the grid search, $D_f$ is chosen to be between 1.8 and 2.8 with a step size of 0.2. $R_m$ varies from 15 nm to 25 nm with an interval of 5 nm, to be both computationally efficient and to avoid the high sensitivity of optical properties to changes in $N_m$. To reach the particle size parameter indicated by the slope of observations at forward scattering, the range of $N_m$ is from 800 to 1500 for $R_m = 15$ nm and from 300 to 1000 for $R_m = 20$ and 25 nm. $R_m$ is fixed in each group of the grid search. The resulting scattering and extinction intensities are evaluated using Equation (2), but only for results in the forward scattering range. The combination of the fitting parameters that provides the smallest $\chi^2$ is defined as the best fit under each scenario with a fixed $R_m$, whose results are given in Table 4 (Iteration 1), with

corresponding scattering results shown with red lines in Figure 8. The three scenarios with $R_m = 15, 20$, and 25 nm generally agree well with observations, and with good agreements in the forward scattering. However, all of them underestimate the backscattering by a factor of three.

In iteration 2 of the grid search, we add in the ice coating. As adding the ice coating onto the aggregates with the same $N_m$ (shown in Table 4, Iteration 1) will change the aggregate size ($R_a$), which strongly affects the forward scattering, $N_m$ is adjusted accordingly. Under each of the three scenarios above with different $R_m$, the particle size parameter is fixed by replacing part of the aggregate material with organic ice coating to reduce the impact on forward scattering. Two parameters ($V_f$ and RI<sub>coat</sub>) are introduced to describe the ice component, and the aggregate radius ($R_a$) is kept constant. $V_f$ changes from 0.7 to 0.9 with an interval of 0.1, and RI<sub>coat</sub> varies from 1.3 to 1.6 with an interval of 0.1, which is summarized in Table 3. Then, the grid search is applied to these two ice coating parameters, and the best fits are found in this iteration. To illustrate the improvement of the fit for backward scattering, we calculate $\chi^2_{\text{back}}$ (value of $\chi^2$ when only backward scattering is considered) and add it to the list of parameters in Table 4.

Using the same criterion of the minimal $\chi^2$ of Equation (2), the best-fit parameter combination is searched under each $R_m$ size scenario, except that under the scenario of $R_m = 25$ nm the value of $\chi^2$ consistently decreases with RI<sub>coat</sub> becoming
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Figure 2, but for the ratio of scattering intensity to UV extinction at varying volume fraction of the aggregate ($V_\phi$) with a spherical coating for a fixed particle size parameter.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Ratio of (a) LORRI to Alice and (b) LEISA to Alice observations for the results from the four iterations listed in Table 4. The red, blue, and green lines indicate iteration 1, iteration 2, and the best fit, respectively, with different line types denoting different monomer radii. Note that the results for iteration 2 and the best fit (green dashed line) are the same when the monomer radius is 20 nm. Altitude-averaged observations are shown as gray error bars. The observation uncertainties include photon noise and read noise (for LORRI and LEISA), which are computed through error propagation when binning the observations to the altitude bins, and also the standard deviation of the observed values in each bin.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Ratio of (a) LORRI to Alice and (b) LEISA to Alice observations for the results of best fits compared with the monodispersed aggregates and bimodal distributed particles from Fan et al. (2022). The gray error bars are the altitude-averaged observations between 7.5 and 47.5 km. The blue asterisks and blue crosses denote results for the bimodally distributed particles and monodispersed aggregates from Fan et al. (2022), respectively. The red boxes emphasize the better agreement of the best fit compared to monodispersed aggregates in backscattering. The bimodal distribution considers two populations of particles. One population is ~1 µm aggregates with a number density of 0.3 cm⁻³ and consisting of ~20 nm monomers; the other is ~80 nm spherical particles with a number density of 10 cm⁻³.</figcaption></figure>

<table><caption>**Table 3**
Testing Ranges of Parameters for Coated Aggregates</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>V<sub>f</sub></td><td>0.7, 0.8, 0.9</td></tr><tr><td>RI<sub>coat</sub></td><td>1.3, 1.4, 1.5, 1.6</td></tr></tbody></table>

smaller to 1.0, with a minimum of 15.90 in the tested range of RI<sub>coat</sub>. There is no local minimum within the physically reasonable parameter range; therefore, we do not consider this scenario with R<sub>m</sub> = 25 nm in the following grid search. In iteration 3, we then fix the values of V<sub>f</sub> and RI<sub>coat</sub> to the results of iteration 2 and further vary N<sub>m</sub> and D<sub>f</sub> with smaller step sizes of 40 and 0.1, respectively, and then conduct iteration 4 with V<sub>f</sub> and RI<sub>coat</sub> changing with a step size of 0.05. After four iterations of such a grid search, we obtain the final sets of parameters that provide the best fit, which are listed in Table 4 with the corresponding χ². The results of the scattering properties of the best-fit haze particles are shown in Figures 8 and 9 in the forms of ratios of LORRI to Alice and LEISA to Alice observations.

After four iterations of the grid search with ice coating added, noticeable improvements can be identified (Figure 8). The set of parameters that corresponds to the best fit is R<sub>m</sub> = 15 nm, N<sub>m</sub> = 880, D<sub>f</sub> = 2.0, V<sub>f</sub> = 0.8, RI<sub>coat</sub> = 1.5. The inclusion of ice coating improves the agreement between model and altitude-averaged observations, especially for backward scattering. However, there are still some problems on the forward scattering. For example, the simulated forward results at the scattering angle of 31.7° are about one time larger than the observed values, and the slope of the IR scattering intensity is slightly larger than that of the observations (Figure 8). In Figure 9, best fits (green lines) are compared with results for monodispersed aggregates (blue crosses) and bimodally

<table><caption>**Table 4**
Values of Parameters Providing Best Fit for Different R<sub>m</sub> and Corresponding χ²</caption><thead><tr><th>Parameter</th><th>Iteration 1</th><th>Iteration 2</th><th>Iteration 3</th><th>Iteration 4 (Best Fit)</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center; font-weight:bold;">R<sub>m</sub> = 15 nm</td></tr><tr><td>N<sub>m</sub></td><td>1100 (100)</td><td>900</td><td>880 (40)</td><td>880</td></tr><tr><td>D<sub>f</sub></td><td>2.0 (0.2)</td><td>2.0</td><td>2.0 (0.1)</td><td>2.0</td></tr><tr><td>V<sub>f</sub></td><td>...</td><td>0.8 (0.1)</td><td>0.8</td><td>0.8 (0.05)</td></tr><tr><td>RI<sub>coat</sub></td><td>...</td><td>1.5 (0.1)</td><td>1.5</td><td>1.5 (0.05)</td></tr><tr><td>χ<sup>2</sup></td><td>1.66<sup>a</sup></td><td>13.26</td><td>12.96</td><td>12.96</td></tr><tr><td>χ<sup>2</sup><sub>back</sub></td><td>...</td><td>9.81</td><td>9.54</td><td>9.54</td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">R<sub>m</sub> = 20 nm</td></tr><tr><td>N<sub>m</sub></td><td>700 (100)</td><td>600</td><td>600 (40)</td><td>600</td></tr><tr><td>D<sub>f</sub></td><td>1.8 (0.2)</td><td>1.8</td><td>1.8 (0.1)</td><td>1.8</td></tr><tr><td>V<sub>f</sub></td><td>...</td><td>0.9 (0.1)</td><td>0.9</td><td>0.9 (0.05)</td></tr><tr><td>RI<sub>coat</sub></td><td>...</td><td>1.3 (0.1)</td><td>1.3</td><td>1.3 (0.05)</td></tr><tr><td>χ<sup>2</sup></td><td>1.40<sup>a</sup></td><td>13.01</td><td>13.01</td><td>13.01</td></tr><tr><td>χ<sup>2</sup><sub>back</sub></td><td>...</td><td>9.96</td><td>9.90</td><td>9.90</td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">R<sub>m</sub> = 25 nm</td></tr><tr><td>N<sub>m</sub></td><td>400 (100)</td><td>...</td><td>...</td><td>...</td></tr><tr><td>D<sub>f</sub></td><td>1.8 (0.2)</td><td>...</td><td>...</td><td>...</td></tr><tr><td>χ<sup>2</sup></td><td>2.82<sup>a</sup></td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Notes. The step sizes of the grid search are shown in parentheses.

a In iteration 1, the calculation of χ² only includes forward scattering.

distributed particles (blue asterisks) from Fan et al. (2022). This suggests that the agreement for ice-coated aggregates is better than that for monodispersed aggregates (χ² = 14.23), especially in backscattering (shown in red boxes in Figure 9(a)), but worse than the fit for bimodally distributed particles (χ² = 12.45). Compared to observations, results from DDA tend to have a smaller backward scattering and larger forward scattering at the scattering angle of 0°–30°. This indicates that the population of smaller-sized haze particles likely exists in the region of ice condensation.
## 5. Discussions and Conclusions

In this work, we introduce inhomogeneous particles to understand Pluto's haze observations. The DDA method is used to obtain the optical properties of coated aggregates, and the multiwavelength observations from the New Horizons spacecraft are used to constrain particle properties. With those done, the grid search approach is employed to find the set of parameters that can best explain the observations. By considering particles as fractal aggregates, a set of parameters, such as $N_m$, $D_f$, and $R_m$, are varied within ranges prescribed based on current knowledge, to evaluate their influence on aggregate optical properties. Most importantly, we consider the influence of organic ice coating materials given the relatively low temperatures detected and the density reduction of organic gases and compute the optical properties of coated aggregates with different coating fractions and refractive indices. Finally, the set of parameters that can best explain the observations is retrieved.

The organic ice coating is a potential candidate to explain observations from the New Horizons spacecraft, and we find that the haze particle having a volume equivalent radius of 0.15 µm and consisting of monomers with a radius of 15 nm results in optical properties that best fit the observations. The volume fraction of organic ice is expected to be ~20% with a refractive index of 1.5. However, such particles cannot explain all the observations at the same time, e.g., forward scattering in the visible. The reason that we achieve better agreement for backward scattering than forward scattering is that the change of the asymmetry factor of the phase function curve caused by the addition of organic ice coating tends to increase the intensity of backward scattering for the particles with the same size parameter. Previous studies show that monodispersed and lognormal populations cannot explain the observed forward and backward scattering at the same time and that a bimodal distribution is needed (Kutsop et al. 2021; Fan et al. 2022), which requires more evidence for the haze production mechanism from observations. In addition, considering the limited number of scattering angles that were observed and that the change of the organic coating refractive index with wavelength is not known, more observational data and laboratory experiments are needed to further constrain the morphology and composition of haze in Pluto's atmosphere.

Pluto's haze particles may be more complex than the ones in our model. A monodispersed population is assumed in this work, but particles are unlikely to be all of the same size in reality. Therefore, distributions of particle sizes, such as a lognormal or power-law distribution, should be considered (Kutsop et al. 2021; Fan et al. 2022). In addition, the centers of the aggregate and the coating are assumed to coincide, and the coating is assumed to be a sphere in this study, which might not hold true in reality since ice can condense on each monomer (Teng et al. 2019). The refractive index of the coating is set as a constant because of the limited knowledge about its composition, but it should change with wavelength similar to that for the aggregates; again, laboratory measurements are required. For the grid search, the choice of the step size in the search for best parameters can affect the results. The range of parameters in our tests was limited owing to the computationally expensive

nature of the simulations. Finally, the presence of coating is deduced from the unexpected low temperature of Pluto's lower atmosphere, but the way the coating is formed remains unclear. Therefore, more observational and laboratory data are vital for an accurate retrieval of haze particles in Pluto's atmosphere.

The authors thank Drs. Maxim A. Yurkin and Alfons G. Hoekstra for their ADDA code. The simulations are conducted in the High-Performance Computing Center of Nanjing University of Information Science & Technology. This work is supported by the National Natural Science Foundation of China (grant No. 42122038). Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (grant No. 80NM0018D0004). Y.L. Y. was supported in part by a NASA NFDAP grant (grant No. 80NSSC19K0823) to Caltech.

## ORCID iDs

Jingyu Wang https://orcid.org/0000-0001-5128-585X

Siteng Fan https://orcid.org/0000-0002-3041-4680

Chao Liu https://orcid.org/0000-0001-7049-493X

Vijay Natraj https://orcid.org/0000-0003-3154-9429

Leslie A. Young https://orcid.org/0000-0002-7547-3967

Yuk L. Yung https://orcid.org/0000-0002-4263-2562

## References

Bi, L., Yang, P., Kattawar, G. W., & Mishchenko, M. I. 2013, JQSRT, 116, 169

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

Dong, J., Zhao, J. M., & Liu, L. H. 2015, JQSRT, 165, 43

Draine, B. T. 2003, ApJ, 598, 1017

Fan, S., Gao, P., Zhang, X., et al. 2022, NatCo, 13, 240

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

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

Gmachowski, L. 2002, ColSu, 211, 197

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

Haynes, W. M. 2013, The CRC Handbook of Chemistry and Physics (94th ed; Boca Raton, FL: CRC Press)

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

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

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

Khare, B. N., Thompson, W. R., Sagan, C., et al. 1989, BAAS, 21, 982

Kutsop, N. W., Hayes, A. G., Buratti, B. J., et al. 2021, PSJ, 2, 91

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

Li, J., Liu, C., Yin, Y., & Kumar, K. R. 2016, JGRD, 121, 3506

Liu, C., Xu, X., Yin, Y., Schnaiter, M., & Yung, Y. L. 2019, JQSRT, 222-223, 170

Liu, F., Yon, J., & Bescond, Y. 2015, JQSRT, 172, 134

Liu, Q. 1997, M OTL, 15, 158

Mackowski, D. W., & Mishchenko, M. I. 1996, JOSAA, 13, 2266

Purcell, E. M., & Pennypacker, C. R. 1973, ApJ, 186, 705

Rannou, P., McKay, C. P., Botet, R., & Cabane, M. 1999, P&SS, 47, 385

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

Teng, S., Liu, C., Schnaiter, M., et al. 2019, ACP, 19, 2917

Tomasko, M. G., Doose, L., Engel, S., et al. 2008, P&SS, 56, 669

Waterman, P. C. 1965, in Proc. of the IEEE, Vol. 53 (Piscataway, NJ: IEEE), 1

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

Wu, Y., Cheng, T., Zheng, L., & Chen, H. 2016, JQSRT, 179, 139

Yee, K. S. 1966, ITAP, 14, 302

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

Yung, Y. L. 1978, ApOpt, 17, 3707

Yurkin, M. A., & Hoekstra, A. G. 2007, JQSRT, 106, 558