Pluto’s Haze Abundance and Size Distribution from Limb Scatter Observations by MVIC

The New Horizons spacecraft observed Pluto and Charon at solar-phase angles between 16° and 169°. In this work, we use the Multispectral Visible Imaging Camera (MVIC) observations to construct multiwavelength phase curves of Pluto’s atmosphere, using the limb scatter technique. Observational artifacts and biases were removed using Charon as a representative airless body. The size and distribution of the haze particles were constrained using a Titan fractal aggregate phase function. We find that monodispersed and log-normal populations cannot simultaneously describe the observed steep forward scattering, indicative of wavelength-scale particles, and the non-negligible backscattering indicative of particles much smaller than the wavelength. Instead, we find it necessary to use bimodal or power-law distributions, especially below ∼200 km, to properly describe the MVIC observations. Above 200 km, where the atmosphere is isotropically scattering, a monodisperse, log-normal, or a bimodal/power law approximating a monodispersed population is able to fit the phase curves well. As compared to the results of previously published articles, we find that Pluto’s atmosphere must contain haze particle number densities an order of magnitude greater for small (∼10 nm) and large (∼1 μm) radii, and relatively fewer intermediate sizes (∼100 nm). These conclusions support a lower aggregate aerosol growth rate than that found by Gao et al., indicating a higher charge-to-radius ratio, upwards of 60e − μm−1. In order to generate large particles with a lower growth rate, the atmosphere must also have a lower sedimentation velocity (<∼0.01 m s−1 at 200 km), which is possible with a fractal dimension of less than 2.


Introduction
Pluto has long been known to have a seasonally varying atmosphere (Elliot et al. 2003;Sicardy et al. 2003). On 2015 July 14, the New Horizons spacecraft flew within 12,500 km of Pluto's surface, revealing the extent of its unique haze. (Stern et al. 2015;Gladstone et al. 2016). Pluto's haze is optically thin, globally distributed, and extends to altitudes of >200 km, with a scale height between 30 and 50 km (Gladstone et al. 2016). Below 80 km, up to 20 distinct layers of between 1 and 10 km thickness are observed over vast expanses of the globe (Cheng et al. 2017). Unfortunately, observational constraints force us to bin our data into 20 km wide altitude bins, such that individual haze layers are not resolved in our analysis. Pluto has been mapped using the Hubble Space Telescope, but only at a single phase angle (<1°), and with a spatial resolution of 600 km (Stern et al. 1997;Buie et al. 2010), whereas the New Horizons flyby has provided the first observations of Pluto over a range of phase angles (15°-170°), and at spatial resolutions of less than <1 km/pixel. In this work we use images acquired by New Horizon's Multispectral Visual and Infrared Camera (MVIC; Reuter et al. 2008) to generate calibrated phase curves of Pluto's atmosphere. We then use these phase curves to constrain the size and distribution of Pluto's haze particles, using a fractal aggregate phase function.
High phase angle images indicate that the atmosphere is strongly forward scattering. This indicates particles with a size parameter of x 1, where x = 2πa/λ, a is the effective radius of the haze particle, and λ is the wavelength of the incident light (for visible light in the range of 400-700 nm, this would necessitate a radius of, a > ∼100 nm). MVIC color composite images, however, reveal that the atmosphere is also very blue ( » 0.5 would require a radius of, a ≈ 10 nm) (Gladstone & Young 2019). This suggests that the particles are fractal aggregates, similar to those formed on Titan (Gladstone et al. 2016). UV solar occultation observations by the Alice instrument (Young et al. 2018) show that the gaseous part of Pluto's atmosphere bears some similarities, in terms of composition, to those of Titan and Triton, with N 2 , CH 4 , C 2 H 2 , C 2 H 4 , and C 2 H 6 all found to be present. In addition, ground-based observations acquired by the Atacama Large Millimeter/submillimeter Array confirm the presence of CO and HCN (Lara et al. 1997;Summers et al. 1997;Krasnopolsky & Cruikshank 1999;Wong et al. 2015;Lellouch et al. 2017). Far and extreme ultraviolet photons drive photodissociation and photoionization of these constituents. With regard to Titan, photolysis produces monomers and other precursors, which then combine and grow as they descend through the atmosphere (Lavvas et al. 2013). The same process is expected to happen on Pluto, albeit with a different charged particle environment, due to the lack of a thermosphere, exposure to the bare solar wind, and the lack of a strong magnetospheric radiation component from a nearby giant planet. Previous works, based on photolysis models, have suggested that aggregates in Pluto's atmosphere have a fractal dimension, D f , of 1.5-2.75, and are composed of 10-20 nm monomers (Lavvas et al. 2013;Gao et al. 2017). Monomer abundance is bounded by the photolysis rates of CH 4 and N 2 , which represent the primary limiting factors in haze production. The vertical haze scattering optical depth is currently only ∼0.004 (Bertrand & Forget 2017). Assuming an average radius of ∼0.2 μm, the haze density near Pluto's surface is ∼0.8 particles cm −3 , or a column mass of 8 × 10 −8 g cm −2 . We use the limb scatter technique, employed extensively for remote sensing of Earth's atmosphere (Cunnold et al. 1973;Flittner et al. 2000;McPeters et al. 2000;Murtagh et al. 2002;Rault 2005;Degenstein et al. 2009;McLinden et al. 2012;Normand et al. 2013) to perform the radiative transfer necessary to model the phase curves of Pluto's atmosphere, as the more familiar bidirectional reflection equations (Hapke 1981;Liou 2002) are inappropriate for our viewing geometries. As shown in Figure 1(a), limb scattering occurs when sunlight enters Pluto's atmosphere, undergoes one or more scattering events between the atmosphere and the surface, then enters the line of sight (LOS) vector between the MVIC and Pluto (McLinden et al. 2012). The tangent altitude is the height measured at the point along the LOS vector which is closest to the surface. The line between the tangent point and the subtangent point on the surface is parallel to the surface normal. Scattering into the LOS vector occurs throughout the atmosphere. For an optically thin atmosphere (τ = 1, where τ is the optical depth), scattering into the LOS vector is dominated by the single scattering contribution from the tangent point (Normand et al. 2013).
The phase curve describes the brightness of a scattering target subject relative to a known illumination geometry, denoted by the phase angle, θ. The phase angle is the angle formed between the incident vector and the emission vector. The incident vector starts at the illumination source, in our case the Sun, and ends at the target, i.e., Pluto's haze, while the emission vector starts at the observer, MVIC, and New Horizons, and ends at the target. The phase angle, θ, is defined as where i is the incidence angle, formed between the incident vector and the normal vector of the target, e is the emission angle, and ψ is the azimuthal angle, which is the angle between the projection of the incidence and emission vectors onto the horizontal plane, defined by the normal vector. Depending on the field, the scattering angle α is sometimes used in place of θ, where α = 180°-θ. In this work, we use the phase angle, θ. The scattering phase function is the first element of the scattering phase matrix, P 11 . In order to be consistent with other works utilizing multiple components of the phase matrix, we will use the notation, P 11 to describe our phase functions. The phase function is normalized, such that The phase function represents the probability that the radiation propagating in a given direction is scattered into an elementary solid angle, making an angle θ with the incident direction (Sharma 2018). It is typical to assume that a particle scatters symmetrically in the azimuthal direction (Sharma 2018), leaving the phase angle as the sole descriptor for the directionality of the scattering event.

MVIC Data
The MVIC uses one of the two feeds of the Ralph telescope, a 75 mm aperture f/8.7 system, with an effective focal length of ∼658 mm (Reuter et al. 2008). The MVIC itself consists of seven independent Charged-Couple Device (CCD) arrays, bonded to a single substrate. Six of the 5024 × 32 pixel CCD arrays are operated in time delay integration (TDI) mode. In this mode, the parallel transfer rate of each row is synced to the relative motion of the image across the detector's surface. TDI increases the effective integration time by 32X, substantially increasing the image's signal-to-noise ratio (S/N). Of the six TDI arrays, there are two panchromatic/unfiltered (MP1, MP2) arrays for broadband high S/N imaging, and four arrays with color filters in the red, blue, near-infrared, and methane windows (MC0, MC1, MC2, MC3). The seventh CCD is a panchromatic frame transfer array, measuring 5024 × 128 pixels (MPF), primarily provided for optical navigation. In all seven frames, the first and last 12 pixels in each row are not optically active. The raw and calibrated data are available via the PDS (Stern et al. 2007). A detailed explanation of our data processing is provided in Appendix A.
We use data only from the Red, Blue, and NIR filters, which have pivot wavelengths of 620 nm, 475 nm, and 878 nm, respectively. We ignore the panchromatic filters, as we are investigating the spectral characteristics of Pluto's haze. We ignore the CH 4 filter (pivot wavelength of 885 nm) because it is a narrow wavelength filter with low S/N, which is encompassed by the wavelength range of the broader NIR filter. Furthermore, we only use images taken with an average pixel scale across Pluto's disk of less than 15 km/pixel. This produces the data set in Table 1, acquired within ∼30 hr of closest approach. Images outside this range cover a redundant range of phase angles but have pixels that are too large to permit probing distinct altitudes in Pluto's atmosphere.
These lower-resolution images are also more susceptible to "off-disk glow". This off-disk glow is an instrument effect, similar to the stray light present in the Long Range Reconnaissance Imager (LORRI) (Cheng et al. 2017). This effect can be seen more than 100 pixels away from the disk, equating to more than 500 km in some images ( Figure 2). To account for this glow, we follow a similar strategy to that of Cheng et al. (2017), whereby we use the airless (i.e., presumably hazeless) body of Charon (Stern et al. 2017) to model and remove the effect. Figure 3 highlights our fit to the data, and provides a comparison with Pluto's off-disk profile. Our technique is detailed in Appendix B.
We present the uncertainty in I/F as the bounds occurring at the 15th and 85th percentile of the pixels in each bin (Appendix B, Figure 18). The 15th and 85th percentile range is equivalent to the bounds at one standard deviation, σ, greater or less than the mean or median for a normal distribution. The variance in our data in a single bin is typically not normally distributed, such that standard deviation is not appropriate for describing out variance. This variability in our binned data, i.e., ∼25%, and filter-dependent (see Appendix B,Figures 17 and 18), is substantially greater than the 3% propagated error described in Appendix C. As such, we do not consider the propagated error, but present our determination of the propagated error above for the sake of completeness in Appendix C.

Phase Curves
We have produced phase curves for Pluto from 0 to 500 km, in 20 km wide altitude bins, at eight solar-phase angles of 16°, Note. Each of the images were captured in the 4 color CCDs. As each CCD is stored in an independent file, each file has a preface relating to the file extensions in Table 1, i.e., XXX_sci.fit would be mc0_XXX_sci.fit, mc1_XXX_sci.fit, etc. All data, i.e., phase, resolution, etc., are given for Pluto, except for those marked with an asterisk, which are given for Charon.  Figure 4 shows an example of our phase curve at 20-40 km in all three wavelengths. Figure 5 shows the relative changes in the high-phase and low-phase components of the phase curves as a function of altitude. The phase curves in the lower atmosphere are strongly forward scattering in all three filters. At 0-20 km, the forwardscattered component is an order of magnitude higher than the back-scattered component in all three filters. A convenient parameter for describing phase curves and phase functions is the asymmetry parameter, often denoted as g (Henyey & Greenstein 1941); g is defined as the first moment of the phase function and quantizes the degree of anisotropy in the scattering from −1 to 1, where g = 1 is fully forward scattering (typified by Mie scattering), g = 0 is isotropic (typified by Rayleigh scattering), and g = −1 is fully backscattering (typified by geometric scattering). We approximate g by fitting a simple monodisperse phase function (3.4.1 Monodisperse and Log-Normal) to the phase curves, and finding the first moment (Liou 2002, Equation 3.4.9(b)). From 0 to 40 km, we find the atmosphere to have a range of g from ∼.6 in the NIR filter to ∼0.7 in the Blue filter. Regardless of scatterer type or distribution, this suggests that the haze is characterized by particles scattering in the Mie regime, where the size parameter x 1. Given the wavelength range of the three MVIC filters, the haze in the lower atmosphere must be characterized by particles on the order of 100 nm, consistent with other forwardscattering measurements (Gladstone & Young 2019;Gao et al. 2017).
As altitude increases, the I/F at all phase angles for the three filters decreases, and they become less forward scattering ( Figure 5). From 0 to 60 km, g decreases slowly, while the total intensity decreases mostly uniformly across the eight phase angles, suggesting that the characteristic particle size is not changing very much, while the total abundance decreases uniformly. This implies that, from 0 to 60 km, aerosol growth has largely ceased, preferential destruction or fallout of particular particles is not occurring, and all scatterer populations are uniformly increasing as they descend to the surface. From 60 to 200 km, g decreases rapidly, as the intensity at . MVIC Red limb scatter profiles of Charon, taken from sequence 0299147977, and of Pluto, taken from sequence 0299127622. These observations were chosen such that they have nearly equal phase angles of, 16°.82, and 16°.15, and the targets have similar angular diameters of 0°.15, and 0°. 19, respectively. The limb profiles were chosen such that the incidence angle of the subtangent points were nearly identical, at 79°.80, and 79°. 85, and the weighted average limb brightness values were nearly identical, at I/F = 0.1513 and 0.151 5, respectively. The fit, given by the red dotted line, was generated for Pluto's limb scatter profile using Equation (6), and the values for MVIC Red in Table 2.  . Approximate high-phase angle (∼165°-170°) and low-phase angle (∼16°-18°) components of MVIC Red, Blue, and NIR phase curves for each altitude bin. θ = 40°, 41°, 166°, 167°, and 169°decreases much faster than the intensity at θ = 16°, 17°, and 18°. This is consistent with the characteristic particle size decreasing rapidly from 60 to 200 km. The difference in the rate of change for the scattering anisotropy from 0 to 60 km, and from 60 to 200 km is shown in Figure 5, in terms of the linear decrease in the low-phase angle component, and the roughly exponential decrease in the highphase component.
The phase curves for the three filters all approach isotropic scattering in the upper atmosphere. The NIR and Red filter turn isotropic at 150 km, while the Blue filter becomes isotropic just below 200 km. This is consistent with the characteristic particle size entering the Rayleigh scattering regime for the NIR and Red filters at 150 km, while the particles are large enough compared to the wavelength range for the Blue filter (475 nm) that they remain in the Mie scattering regime until 200 km. If we take the transition from Rayleigh to Mie Scattering to be X = 0.2, (Wiacek et al. 2013), this would suggest that the characteristic particle size is ∼30 nm at 150 km, and ∼15 nm at 200 km. This is about as small as the fractal aggregate haze particles can be, if they are composed of monomers with spherical radii of 5-10 nm . In light of this, and the constant mild decrease in I/F for the phase curves in all three filters, it is likely that minimal aerosol growth occurs between 200 and 500 km, while the total abundance of particles grows steadily.

Radiative Transfer
A full radiative transfer description of scattering from Pluto's limb would include ray tracing through spherical shells, and would account for multiple scattering (Bourassa et al. 2008). These multiple scattering events are illustrated in Figure 1(b). Since none of our observations intersect the surface of Pluto, we ignore all scattering between the surface and the MVIC (Sun → surface → MVIC, Figure 1(b) A and Sun → atmosphere → surface → MVIC, Figure 1(b) D). The results of Fan et al. (in preparation) indicate that the upper limit of the contribution from secondary scattering initiated by the ground (Sun → surface → atmosphere → MVIC, Figure 1(b) E) accounts for less than 10% of the observed radiance at the wavelength of the three MVIC filters used in our work ( Figure 4). For this reason, we ignore secondary scattering initiated by the ground, together with higher-order multiple scattering. We also ignore secondary (and multiple) scattering initiated by the atmosphere, (Sun → atmosphere → atmosphere → MVIC, Figure 1 We determined the magnitude of the scattering due to Pluto's primary gaseous components, N 2 and CH 4 , using the particle number density described in Young et al. (2018), and the Rayleigh scattering phase function. At 370 nm, where laboratory data for N 2 and CH 4 cross-sections exists, we find that the Rayleigh scattering at the surface (0 km) due to N 2 is I/F ≈ 10 −3 , while Rayleigh scattering due to CH 4 is I/F ≈ 10 −7 . Given that Rayleigh scattering intensity approximates to I ∝ 1/λ 4 , the intensity at 370 nm is three times stronger than for our Blue filter pivot wavelength of 492 nm. This also means that Rayleigh scattering is even further reduced in the Red and NIR filters.
As described in detail in Appendix D, we ignore all scattering from Pluto except for single scattering from Pluto's haze, which gives us

Grid Search and Statistics
Our goal is to find a phase function that matches the observed limb scatter phase curves at all three wavelengths simultaneously at each altitude bin, for all altitude bins. Between these, we are looking to match the phase angle dependence, relative brightness difference between high and low phase, and slope d(I/F(λ))/dθ. To do this we used an R 2 grid search, where each element of the grid is Here, SSE is the sum of squared error (sometimes called the sum of squares of residuals) and SST is the total sum of squares, each given by where C is the total number of phase angles in our phase curves, and k is each point in both the observed phase curve I/F obs (θ), and the phase function, P 11 (θ). Our objective is the distribution and abundance of the scattering population, which are the driving aspects of the shape of the modeled phase curve. Given that ( ) ( ) s l N s s ,

T T
are independent of phase angle, under the assumption that the aerosols are randomly oriented, we replace them in Equation (5) with a scaling factor, defined as R 2 ideally has a range of 0-1, with a value closer to 1, indicating that a greater proportion of variance is accounted for by the model, and that the model more accurately predicts the data. Here, R 2 0 indicates that the mean of the data, i.e., a horizontal line, provides a better fit to the data than the model. Negative R 2 values are possible, but typically occur for equations that do not contain a constant term.
We plot our data using the median, and report the variation for the 15th and 85th percentile values. However, we conduct our fitting routines using a Monte Carlo (MC) simulation, with at least 320 runs, to draw our phase curve values for each wavelength from the distribution of pixels in each of the altitude and phase bins. For each MC run, we redo our grid search, and locate the best-fit parameters. At the end, we determine our best-fit parameters by taking the average of the MC runs, and report our uncertainty with the 15th and 85th percentile values for each fit parameter.

Size Distribution
Background information about the phase functions and the choices we made is provided in Appendix E. Phase functions produced by a distribution of particle sizes, characterized by the differential particle radius number density distribution, n(r), are found by integrating the phase function over all particle sizes (Grainger 2020): where N p (r) is the number of particles with radii between r and r + dr, σ s and σ e are the scattering and extinction cross-section of the particle, and r is the effective particle radius; a UV is the normalization constant necessary to find the total number of particles N 0 , and to match the extinction coefficient observations made by the New Horizons Alice UV spectrometer, determined by is the extinction coefficient of the haze observed by Alice during the solar occultation of Pluto at the tangent altitude of s T (Young et al. 2018, Appendix B). In order to do this, we needed to extrapolate the complex index of refraction from Ramirez et al. (2002), from 200 nm to 165 nm (the wavelength of the Alice occultation observations); we did so linearly. Note that from here on, P 11 will refer to the phase function produced by a distribution of scatters, while P 11,Individual will refer to the phase function produced by a single scatter of characteristic size between r and r + dr.
The total number of particles per unit volume, N 0 , is given by for spheres, and for fractal aggregates. Here, Q s,e denotes the scattering and extinction efficiencies, respectively. The effective particle radius for a spherical particle is simply the radius of the sphere, whereas for a fractal aggregate, the effective radius, R f , is given by where r m is the radius of a monomer, N m is the number of monomers in a fractal aggregate, and D f is the fractal dimension. For fractal aggregate distributions, we assumed a particle size range between 2 monomers (R min = r m 2 1 /D f ) to R max = 1 μm. The minimum size ensures that the particle is actually an aggregate, and not a single monomer. The maximum size was initially chosen based on the upper range used by Gao et al. (2017). Testing of larger R max , up to 1 mm, revealed that the particular slope of the high-phase component of the phase curves necessitates that there be a dramatic lack of particles larger than 1 μm ( Figure 8). This will be discussed in greater detail in 3.4., Population Types.
We also assume that the fractal dimension is D f = 2. This is because only the fractal dimension of D f = 2, has been rigorously tested, by Tomasko et al. (2008). Unpublished results from further testing suggest that the code can produce accurate phase function for 1.5 < D f < 2.5. However, these may also produce physically irrelevant albedos, and/or polarization values (Email communications with M. Lemmon 2020). As will be discussed in 3.4. Population Types, we investigated this range of D f as the fractal dimension, as an important parameter, controlling the fallout time of the haze.
We assume that monomers maintain a constant size after formation, and through fractal aggregation; as such, we fix r m as a constant for all altitude bins. The final P 11 must fit the phase curves of all three MVIC filters. Throughout our tests, we found that modeled phase functions with different r m values matched the observed phase curves equally well for all r m 10 nm, while the quality of the fits diminished as r m increased, for r m > 10 nm. As will be discussed below, we found that the 10 nm monomer fractal aggregate distribution reported by Gao et al. (2017) produced phase functions which matched our phase curve better than the 5 nm monomer fractal aggregate, the 10 nm spherical, or the 5 nm spherical distributions given in Gao et al. (2017). Given this observation, together with our other test, we assumed that r m = 10 nm. In Section 4.1 we will discuss the fractal dimension for different monomer sizes, with specific reference to the fractal dimension, D f .

Population Types
We tested several particle population types, including monodisperse, log-normal, bimodal, and power law; these are each defined in their relevant sections. For each population type, we tested the phase function produced by every pertinent variation of size and distribution of particles (Equation (7), ) against the observed I/F phase (Equation (7), . At each altitude bin, and for each population type, we identified the specific combination of size and distribution required to produce phase functions which best matched the observed phase curves, by determining which had the highest R 2 value (3.2 Grid Search and Statistics), verifying each via graphic interpretation. We then compared the different populations against one another to determine which were relevant to Pluto's atmosphere.
We will refer to three distributions as small, intermediate, and large-size scatterers. This loose distinction is useful for referencing the distributions in terms of which population category is responsible for a particular aspect of the phase curve. While not strictly in the traditional parameter-size cutoff (Wiacek et al. 2013) we define these groups as Rayleigh-size scatterers, ( ) R Int f ≈ 100-500 nm, and Mie-size scatterers, ( ) R Mie f ≈ 500 nm-1 μm.

Monodisperse and Log-normal
A monodispersed population is characterized by a single particle size. A log-normal distribution is given by Where m R f is the mean, and s R f is the standard deviation of the natural logarithm of the effective radius, R f . Gao et al. (2017) determined n(r) as functions of altitude and particle radius for spherical and fractal aggregate particles, composed of monomers with r m values of 5 and 10 nm, respectively, using the 1D community aerosol and radiation model for atmospheres (CARMA). We found that the log-normal distributions produced by Gao et al. (2017) were unable to simultaneously fit the steep forward-scattering component, and the non-negligible back-scatter component present in Pluto's atmosphere from 0 to 200 km ( Figure 7). As expected, a monodisperse population of scatterers was similarly unable to match the shape and slope of the phase curves ( Figure 7). Of the distributions modeled by Gao et al. (2017), we found that the 10 nm fractal aggregate population best matched our phase curves, being able to match the difference in magnitude between high phase and low phase. This is likely to be because below 150 km, the 10 nm distribution has the widest particle size standard deviation, ) R Ray f occur in sufficient abundance to produce phase functions with the right order of magnitude difference between high phase and low phase. Since the 10 nm fractal aggregate fits better than any other distribution, and because our many tests on all the other population types revealed that r m 10 nm produce phase functions with a nearly identical quality of fits, we use r m = 10 nm for all populations going forward.
The problem remains, however, that there is not the correct disparity in abundance between ( ) R Mie f and ( ) R Ray f to produce the appropriate slope at high-phase angles. These distributions fail because the uniquely steep forward-scattering component necessitates large fractal aggregate scatterers. Large ( ) R Mie f , with the appropriate slope, have low-phase angle components with at least two orders of magnitude lower I/F, which is not what we observe ( Figure 4). Conversely, to produce the low-phase component, which at most is only one order of magnitude lower than the high-phase component (0-60 km, Figure 5), requires smaller sized ( ) R Ray f scatterers. These ( ) R Ray f , however, scatter more isotropically, and are unable to produce steep slopes at high phase. To demonstrate this, we found best-fit phase functions for only the high-phase component of our phase curves, θ > 90, and low-phase component of our phase curve, θ < 90. Figure 6 demonstrates the need for a population of diverse particles, by showing the best-fit results for the Blue filter. While not depicted in Figure 6, the same is true for the MVIC Red, and MVIC NIR filters. Figure 6 shows that the best fit, using monodisperse aerosols, is unable to reproduce the observations in the three MVIC filters.
We also investigated other log-normal distribution by modifying the mean effective radius, m R f , and the standard deviation, s R f , of the Gao et al. (2017) distribution. No amount of manipulation was able to produce a size distribution with a phase function to match the shape of our phase curves appreciably better than the unmodified distributions of Gao et al. (2017).

Bimodal
The bimodal distribution is a population characterized by two scatterers, ( ) R Big f , and R f (Small), with arbitrary weightsw(Big) and w(Small)-that summed to one. We utilize the conversions from Grainger (2020) for the integration of discrete particle sizes to a summation of weights. We found that the bimodal population was able to match the shape and slope of the phase curves very well at all altitude bins (Figures 7(a)-(c)). From 0 to 180 km, the ( ) R Big f value is between 500 nm and 1 μm ( is between 20 and 100 nm (∼ ( ) R Ray f ). From 0 to 200 km, the small particles are at least three orders of magnitude more heavily weighted than the large scatterers ( ( ) ( ) ) = = w w Small 0.99, Big 0.01 . Above 200 km, the weights start to deviate randomly, varying from a uniform distribution ( ( 0.999, Big 0.001 . At this altitude range, however, the size of ( ) R Big f and R f (Small) are very close to one another, and both are in the R f (Ray) size range. This indicates that above 200 km, the population is very similar to a monodisperse population, characterized by particles small enough to scatter isotropically at all wavelengths, as shown in Figure 5.
In Tables 3-6 we include the best-fit parameters for the size of the particles, the weights, and the approximation of the total number of particles the weights suggest, given the scattering cross-section of the particles. We approximate N with   Best-fit phase functions are produced for the power-law, bimodal, trimodal, and monodisperse distributions, and compared against the phase functions produced by the distribution given in Gao et al. (2017). Our observations from MVIC Red, Blue, and NIR filters are presented as filled circles, with 15th and 85th percentile error bars representing the variation in the data bins. All populations are generated with r m = 10 nm and D f = 2 for the 140-160 km altitude bin. The dashed lines are the median (i.e., best fit) of the modeled phase functions, while the shaded area around the dashed lines indicates their 15th and 85th percentiles. As the altitude increases, it becomes more difficult for the power law to match the steep high-phase component. The bimodal and trimodal distributions do not share this problem. RIGHT) R 2 values, showing how well the models (lines) fit the data (circles) in the figures to the left, for each altitude bin from 0 to 100 km, where 1 is the best score. Error bars show the 15th and 85th percentiles of the score. As the altitude increases, the bimodal and trimodal distributions appear to fit the data best.  Best-fit phase functions are produced for the power-law, bimodal, trimodal, and monodisperse distributions, and compared against the phase functions produced by the distribution given in Gao et al. (2017). Our observations from MVIC Red, Blue, and NIR filters are presented as filled circles, with 15th and 85th percentile error bars representing the variation in the data bins. All populations are generated with r m = 10 nm and D f = 2 for the 340-360 km altitude bin. The dashed lines are the median (i.e., best fit) of the modeled phase functions, while the shaded area around the dashed lines indicates their 15th and 85th percentiles. All of the modeled phase functions appear to fit the observations equally well. RIGHT) R 2 values, showing the how well the models (lines) fit the data (circles) in the figures to the left, for each altitude bin from 0 to 100 km, where 1 is the best score. Error bars show the 15th and 85th percentiles of the score. The R 2 for all the populations modeled in this study are nearly identical, suggesting that at this altitude, all of the distributions are approximating a monodisperse population.

Power-law
We describe ( ) n R f using a power-law distribution, where and b is the power-law exponent.
As can be seen in Figures 7(a)-(c), the power law matches the shape and slope of the phase curves nearly as well as the bimodal distribution. In fact, a deep dive into the data reveals that, more often than not, the size and abundance of the bimodal distribution will intersect the two points on the power law. The key difference between the two types is the presence of R f (Int) in the power-law distribution, while the bimodal distribution essentially sets this abundance to zero. This indicates that R f (Int) values are not a defining factor of our phase curves.
The power-law exponent is at its lowest at 20-40 km, where b = 3.5. This is the most equally distributed of the power-law populations, and still there are five orders of magnitude more scatterers with size R f (Ray) than R f (Mie). As we move to higher altitudes, b steadily increases, indicating a greater disparity between R f (Ray) and R f (Mie). The power law becomes steepest, b ≈ 6, at 200 km, where the atmosphere transitions to an isotropically scattering form. This is the same altitude at which we saw the bimodal population begin to approximate a monodisperse population. This again suggests that, above 200 km, the atmosphere is dominated almost entirely by small, isotropically scattering particles.
As a side note, the very steep and dramatic change in the power-law exponent at 180 and 200 km, is likely to be the originate from low signal and low S/N observations, which were not removed, as they were not yet negative. These skewed the grid searches to try and fit such implausibly dim observations. Further investigation revealed that power-law exponents of b ≈ 5.5 (the approximate mean of b at 160 and 220 km) produced phase functions with similar R 2 values, and matched the shape of the phase curves equally well upon inspection.

Trimodal
The bimodal and power-law distributions suggest the same story for the relative abundance of R f (Ray) and ( ) R Mie f . They imply a different story for R f (Int), however, where the bimodal imposes a complete absence, and the power law imposes that R f (Int) is at the average of the logarithmic abundances for R f (Ray) and ( ) R Mie f . To try and determine the importance of R f (Int), we investigated a trimodal distribution. Similarly to the bimodal distribution, the trimodal distribution is characterized by three scatterers, Small that summed to one. These weights are converted to physical abundances using Equation (18). Altogether, the trimodal distribution resulted in a six dimensional grid search, with 214305 unique combinations and corresponding phase functions. Unlike the other distributions, we investigated all of them by sorting them using their R 2 value to determine whether a pattern emerged to suggest any principal components or important characteristic. This resulted in a distribution of R 2 , with a discernible mean and deviation. We defined our "best-fit" parameter and phase functions as those with an R 2 above the 99.7th percentile.
The best-fit trimodal parameters all had an R f (Int) abundance consistent with either a bimodal distribution, or a power law.
This means that it is not possible to immediately remove the degeneracy between power-law and bimodal distributions. As can be seen in the difference between the bimodal and powerlaw phase curves in Figure 7, R f (Int) controls the shape of the phase function in the range where we lack color observations (40°< θ < 165°). This does reveal, however, that R f (Int) must have a relative abundance as compared to R f (Ray), i.e., less than the abundance predicted by a power law. As with the power-law and bimodal distributions, we also note that, for all altitude bins, the most abundant populations of scatterers must lie in R f (Ray). This is not found in the monodisperse or lognormal distributions, which have a μ Rf at ( ) R Int f in the lowest altitude bins. Here, they can match the difference between the low-and high-phase angle brightness, but not simultaneously with the high-phase angle slope (Figures 6 and 7).
As can be seen based on the values given in Table 5, the R f (Med) in all altitude bins was very close in value to either R f (Small) or R f (Big); in fact, they represented the next available value, being as close to either R f (Small) or R f (Big) as the grid search would allow. This would seem to indicate that the trimodal distribution is trying to approximate bimodal, rather than power-law distribution. While this is plausible, it should be viewed with skepticism, as we do not have phase coverage such that R f (Med) being distinctly set apart from R f (Small) or R f (Big) would make a difference in the fits.
The trimodal investigation also revealed that while the size and abundance of the scatterer at R f (Ray) was mostly consistent, there was some variability and multiple combination of size or abundance that would fit the phase curves well. This was not the case for the scatterers at R f (Mie), however, where every single one of our best fits had the same size and abundance. For example, at 20-40 km, it was crucial that the trimodal fit contain a scatterer, with R f ≈ 750 nm and w = 0.01 (N = 491,972). Whether this particle was represented by R f (Med) or R f (Big) was not critical, so long as one of them had this particle size at this particular abundance. We saw the same behavior for the rest of the altitude bins. This is likely to be due to the slope at high-phase angles, which is strongly controlled by the size of particles at R f (Mie), while the difference between the high-and low-phase components requires a finely controlled abundance of particles at R f (Mie).
As noted above, occasionally R f (Med) would occupy this crucial particle size at R f (Mie). In these cases, ( ) R Big f always diverged from a power law in size and abundance, so that ( ) R Big f was less abundant than a power law would predict. Of these best fits where R f (Med) occupied the crucial spot, the best of the best saw ( ) R Big f occupying a size and abundance close to R f (Med), indicating more strongly that there is a dramatic decrease in the abundance of scatterers larger than the crucial particle size. This means that while we cannot definitively say whether or not the distributions follow a bimodal or a powerlaw distribution at ( ) R Int f , we can say that for scatterers larger than the scatterer size responsible for producing the slope at high-phase angles, there must be a precipitous decline in abundance. This can be further confirmed by adjusting the particle range tested for power-law distribution. We initially set R fMax = 1 μm to be consistent with the modeling preformed by Gao et al. (2017). By adjusting R fMax for the power-law distribution, which, unlike the bimodal and trimodal varieties, is a continuous distribution, we saw that the phase functions dramatically change in the shape of their slope at high-phase angles and/or the difference in brightness between high-and low-phase angles ( Figure 8). All of this means that it is necessary to have a specific maximum particle size at the right abundance, and then to have an extreme paucity of larger particles, even if this involves a discontinuity in the particle abundance distribution. An in-depth exploration of maximum particle size actually reveals that the best maximum particle size is R fMax = 1.4 μm. We did not use this value in our analysis, however, because the increase in goodness of fit was minor, and R fMax = 1 μm is easier to compare across the other distribution types.

Three Regimes of Populations
The results of our population investigation support the empirical analysis of the high-and low-phase angle brightness performed in 2.2., Phase Curves. From 500 to 200 km, the atmosphere can be characterized by a tightly distributed population of scatterers at R f (Ray). In this altitude range, photolysis is still taking place (Bertrand & Forget 2017).
From 200 to 60 km, the atmosphere is characterized by two dominant populations of scatterers. The first is a population at R f (Ray), several orders of magnitude more abundant than any other population. The second is a population at R f (Mie), which is far less abundant, but still crucial, as it is required to produce the appropriate phase functions. As can be seen in the first few columns of Tables 4 and 5, the size of the fractal aggregates remains roughly the same for this altitude range. Here, the minimum particle size is about 30-50 nm, and the maximum particle size, referred to as the crucial particle in Section 3.4.4 Trimodal, is about 500-750 nm. While the particle sizes remain the same, the abundances for both increase inversely with altitude, in line with the occultation measurements made by Alice (Young et al. 2018). In order for the lower 200 km to be characterized by a bimodal distribution following this pattern, the scatterers at R f (Mie) would need to grow in abundance independently of R f (Int). It is difficult to imagine a mechanism that would preferentially produce only R f (Mie) scatterers at all altitude bins. That being said, the phase curves and R 2 value in Figure 7 clearly show that bimodal and trimodal populations, with their paucity of R f (Int), are better able to fit the phase curves, and especially the steep slope at high phase angles.
As can be seen in Figure 9, the power-law exponent steadily decreases until 60 km, where it begins to turn over, then increases in the lowest altitude bin at 0-20 km. A similar trend is seen in the parameters for the bimodal and trimodal distributions, and indicates that the abundance disparity between R f (Ray) and R f (Mie) steadily decreases from 200 to 60 km; however, starting at 60 km, this disparity stops decreasing, and briefly shows signs of increasing. This is the same altitude range where a dramatic change in temperature was observed by the REX instrument and ground-based stellar occultations (Sicardy et al. 2016;Hinson et al. 2017).
These three regions of aerosol growth (0-60 km, 60-200 km, and 200-500 km) seem to occur at roughly the same altitudes where the number density of acetylene (C 2 H 2 ), ethylene (C 2 H 4 ), ethane (C 2 H 6 ), and hydrogen cyanide (HCN), as well as the temperature, all experience dramatic slope changes (Lellouch et al. 2017;Young et al. 2018). This indicates either that the processes controlling the abundance of these gases also controls the growth rate of the aerosols, or that these gaseous components are directly responsible for aerosol growth via condensation. These regions also coincide with the termination of photolysis at ∼200 km (Bertrand & Forget 2017), and a dramatic shift in the abundance of haze precursor material and haze (Figure 9).

Fractal Dimensions
We have until now assumed a fractal dimension of D f = 2 for all distributions, for the sake of consistency with Tomasko et al. (2008), and Gao et al. (2017). However, D f plays an important role in determining the porosity of the fractal aggregate, which in turn dictates the sedimentation velocity, cross-sectional area, and coagulation rate. For all the distribution types discussed in Section 3.4 , Population Types, we investigated D f = 1.5, 2, 2.5, while keeping all other parameters as originally described.
A lower D f results in lower sedimentation speed (Figure 10), increasing the time particles have to interact with one another and potentially grow, either by condensation or agglomeration. Although the model is not rigorously tested for a fractal dimension other than 2, as noted above, it will accurately show the difference between the monomer phase function and the aggregate. That is, a lower D f will lead to a less compact aggregate, which will show in the forward scattering.
While the model does account for interference effects within the fractal aggregate, it is also partially parameterized, and may therefore fail to accurately predict polarization and albedo values for aggregates with large numbers of monomers, especially when D f > 2. Given that we focus exclusively on the shape of the phase function in this work, these shortcomings are not important. However, work should be done in the future to ensure that any phase functions with D f ≠ 2 that fit the observed phase curves are physically possible (Email communications with M. Lemmon 2020).
Changing the fractal dimension did not affect the quality of the fits or the R 2 value in any appreciable way. For the bimodal and trimodal distributions, the best-fit sizes for R f also remained the same. As can be seen in Figure 11 for the power-law distribution at 20-40 km, however, changing the fractal dimension resulted in a change in the relative abundance disparity between ( ) R Ray f and R f (Mie), and the total abundance ( ) n R f . This is because the fractal dimension directly affects the extinction cross-section, σ e , where a larger D f results in a larger σ e . As we normalize our size distribution equations, using Equation (11), and the Alice UV extinction coefficient, β Young , we require that a decrease in σ e must be compensated for by an increase in the particle abundance, and vice versa. The biggest compensation occurred for the R f (Mie) scatterers, Figure 9. Best-fit power-law exponents as a function of altitude for D f = 2.0, r m = 10 nm (Blue).The upper x-axis shows the asymmetry parameter, g, for a phase function produced at MVIC Blue wavelengths by a fractal aggregate with the aforementioned fractal dimension and monomer size. (Green) The mass density of precursor material produced by photolysis, as predicted by Bertrand & Forget (2017). (Orange) The mass density of haze produced through aerosol growth of photolyzed material, as predicted by Bertrand & Forget (2017).  where, at 20-40 km, we observed an abundance difference of four orders of magnitude between D f = 1.5, and D f = 2.5.
Meanwhile the R f (Ray) abundance only changed by about one order of magnitude over the fractal dimensions. This is likely to be due to a more significant change occurring in σ e for larger scatterers than for smaller ones. For all D f tested, we noticed the same transitions at the same altitude as discussed in Section 4.1, Three Regimes of Population, for all population types.

Growth Mechanics and Implications
The abundance of R f (Ray) observed in the bimodal and power-law distributions suggests that the coagulation rate is lower than that predicted by Gao et al. (2017). This may be due to the haze having a higher charge-to-particle radius ratio, e − /R f . Figure 5 (top) of Gao et al. (2017) shows variations in the particle size distribution 1.6 km above the surface of Pluto, together with changes in the particle charge-to-radius ratio, given in e − μm −1 . They report a charge-to-particle radius ratio of 30e − μm −1 , which produces a distribution of fractal aggregate in agreement with the retrieved mean particle size (∼0.1-0.2 μm) based on the forward-scattering measurements (Gladstone et al. 2016), which is 2-4 times that obtained from observations of Titan aerosols (Lavvas et al. 2013;Larson et al. 2014). Following the pattern shown in their Figure 5, and drawing a correlation to the high R f (Ray) abundances in our power-law and bimodal distributions, we suggest that the haze has an e − /R f 60e − μm −1 . This is further supported by Figure 5 (bottom) of Gao et al. (2017), which shows profiles of extinction coefficients, β e (Gao et al. 2017 uses α), corresponding to the particle charge-to-radius ratios of Figure 5 (top). While this reveals that β e is not significantly perturbed when the charge-to-radius ratio is varied, it is also obvious that a higher e − /R f results in a lower extinction coefficient below 150 km. While the profiles of β e nearly match the observed β e from the Alice UV solar occultations, they overestimate β e below 150 km, even for the largest charge-to-particles radius ratio tested by Gao et al. (2017), e − /R f = 60e − μm −1 . Figure 5 (bottom) suggests that in order to matchβ e below 150 km, e − /R f must be greater than or equal to 60e − μm −1 .
In Gao et al. (2017) Figure 5 (top), population distributions with a high e − /R f have a high abundance of scatterers at R f (Ray), but a complete absence of R f (Mie) scatterers. As we have shown in Figures 7, 11, and 12, it is necessary have an abundance of R f (Mie) scatterers greater than the abundance predicted by Gao et al. (2017) in order to simultaneously match the difference in brightness between high and low phase, the slope at high phase, and β e observed using Alice. This means that while a lower coagulation rate is necessary (possibly through an increased e − /R f ), there must also be a mechanism which allows the aerosol to grow in size to R f (Mie), and to a sufficient abundance to match the shape of β e . One possibility is to decrease the sedimentation velocity of the haze by decreasing the fractal dimension. Using the sedimentation velocity equations from Gao et al. (2017) we show in Figure 10 that aerosols of R f 100 nm and D f = 1.5, descend through the atmosphere an order of magnitude slower than aerosols with D f = 2. This decrease in velocity gives the aerosols more time to coagulate and form larger particles, even with the larger e − /R f .
In order to describe the observed phase curves, ( ) R Ray f must be present throughout the entire atmosphere in high abundance, and there must be R f (Mie) particles orders of magnitude more abundant than determined by Gao et al. (2017) (Figures 11 and  12). Without a more complete phase curve, we are unable to determine the necessary amount of R f (Int), although we can place upper and lower limits on the abundance. We were able to determine that the abundance of ( ) R Int , f must be less than, or about equal to, what would originate from a power-law distribution. Without a mechanism that preferentially destroys or removes ( ) R Int , f we can also assume that the abundance of R f (Int) must be greater than, or about equal to, the abundance of ( ) R Mie f emanating from a bimodal distribution. It is, however, possible that a mechanism does exist to preferentially remove ( ) R Int , f which point to our distribution being bimodal rather than a power law. This will require further modeling outside the scope of this paper.
The distribution and growth pattern of the aerosols could be indicative of three regimes of aerosol growth, controlled respectively by condensation, UV photolysis, and temperature. Figures 5, 10 and 12 all demonstrate areas where the growth rate and primary population type changes; at ∼60 km, and at ∼200 km. These transitions are highlighted in Figure 13, which shows the logarithmic change in abundance as altitude decreases for three particle sizes at the characteristic sizes of ( ) ( ) R R Ray , Int , f f and R f (Mie). From 500 to 200 km, all three particle sizes experience little increase in abundance, with the smallest particles experiencing the greatest abundance increase, and the largest experiencing the least. This indicates that there is little aerosol growth. From 200 to 60 km, the large particles begin to rapidly increase in abundance, while the small particles stop increasing in abundance, reaching stagnation at ∼150 km. This indicates that larger particles are starting to grow from the smaller particles, utilizing the small particles as fast as they can be delivered from above. From 60 to 0 km, all three particle sizes then converge, and increase in abundance at the same rate. This, along with the return to small particles increasing in abundance, suggests that aerosol growth has largely ceased. 200 km in Pluto's atmosphere has been shown to be the altitude of peak condensation of acetylene (C 2 H 2 ), ethylene (C 2 H 4 ), and ethane (C 2 H 6 ) by means of heterogeneous nucleation (  above 400 km, there are too few aerosols for condensation to be important, and below 200 km, the temperature is too high to permit condensation. We point out that Wong et al. (2017) show that condensation occurs from 500 to 200 km, and abruptly terminates at 200 km. A similar trend at these altitudes occurs with the peak and termination of CH 4 photolysis (Bertrand & Forget 2017 Figures 2-4). Bertrand & Forget (2017) attribute CH 4 photolysis to the production of haze precursor material, which steadily increases from 500 km, accelerating as it peaks at 200 km, before completely terminating at just below 200 km. In addition to, and as a consequence of, "C 2 H x " condensation and the UV photolysis of CH 4 to make precursors, the temperature in Pluto's atmosphere experiences several transitions of its own, as observed by REX and Alice (Hinson et al. 2017, Young et al. 2018. These three aspects of Pluto's atmosphere follow the same path as the population distributions. It begins with a mostly monodisperse population, giving rise to isotropic scattering, before rapidly evolving into disparate populations, with an abundance of R f (Ray), and the a maximum particle size of 500-750 nm. On this basis, we suggest that aerosol growth is limited above 200 km, and produces mostly precursor material at R f (Ray). Where this haze nucleation and condensation is occurring, the temperature in Pluto's atmosphere is roughly constant at around 70 K. At 200 km, the spike in precursor production and condensation should lead to the formation of an abundance of aerosols, mostly at R f (Ray), but at other sizes as well. Below 200 km, where condensation and photolysis has ceased, the temperature begins to increase from ∼70 K at ∼300 km to ∼110 K at ∼40 km. In this more kinetically active environment, free from condensation, the haze particles rapidly interact to form larger particles by means of fractal coagulation. The coagulation rate is such that R f (Mie) increase in abundance by about five orders of magnitude between 200 and 50 km. In this same range, the R f (Ray) particles barely increase by one order of magnitude, showing that the coagulation rate of particles is slightly less than the production rate of precursor material (Figure 13). This temperature-controlled coagulation is further supported by the observation in Figure 13 that at 60 km, the coagulation rate starts to diminish, and all populations uniformly increase just before the surface. This is consistent with the temperature inversion observed by REX and Alice at the boundary surface layer (Young et al. 2018). Here, the temperature first approaches an inflection at ∼110 K, before rapidly cooling to ∼40 K at the surface. In this kinetically less active environment, we would expect the aerosols to cease interactions and proceed to fallout.
These hypotheses suggest possible formation mechanisms for the bimodal or power-law distributions that must be present to produce the phase curves observed by the MVIC. However, further microphysical haze and general circulation modeling will need to be conducted to test these assertions. In addition, more complete radiative transfer limb scatter analyses, which could include multiple instruments that are able to fill the phase coverage gap, will be needed to break the degeneracy of a bimodal or power-law distribution. Breaking this degeneracy will favor one formation mechanism over another.

Summary and Conclusions
The main results of this work are: 1. We have characterized and removed an instrumental bias, which added spurious brightness far off the disks of Charon and Pluto. This off-disk glow accounts for up to 30% of the observed brightness in the lower atmosphere, and greatly distorts the shape of the phase curves, particularly in the back-scatter regime. 2. We have produced phase curves of Pluto's atmosphere, based on the MVIC observations, using the limb scatter technique. These observations were corrected for the above described off-disk glow. The phase curves were produced for the MVIC RED (pivot wavelength of 620 μnm), MVIC BLUE (475 nm), and MVIC NIR (878 nm) filters. Our data set covered an altitude range of 0-500 km, as well as low and high solar-phase angles of 16°.15 to 18°. 24, 38°.85, and 165°. 96 to 169°.36. For this work, we binned all the data which matched our subselection criteria, as described in Section 2.2, Data Processing, into 20 km wide altitude bins, and onedegree width phase bins. This produced the data sets given in Appendix B, Figures 17 and 18. 3. The phase curves show that Pluto's lower atmosphere (<200 km) is strongly forward scattering, indicative of a population of scattering particles in the Mie regime, meaning that the size of the particles is approximately equal to the wavelength of light. In the lower atmosphere, there is a non-negligible amount of back-scatter. This indicates that there must be a distribution of particles, as the phase functions produced by only the aforementioned Mie particles have almost no back-scatter component. As such, there must be an additional population of small size scatterers, which occupy the Rayleigh regime, and are much smaller than the wavelength of light. In the upper atmosphere (>200 km), scattering becomes nearly isotropic. This suggests that between the upper atmosphere Figure 13. Log abundance gradient for characteristic particle sizes of the bestfit power-law distribution in Figure 12. We omit the 200-220, and 220-240 altitude bins for R f (Int) and ( ) R Mie f , as their solutions are unconstrained. The yellow, orange, and purple regions denote the three regimes of aerosol growth. Yellow: 500-200 km. Here, photolysis is active, and little aerosol growth occurs. Orange: photolysis has ceased, and aerosols grow rapidly, and increase in abundance. Purple: the temperature drops, and aerosol growth decreases precipitously, all sizes increase in abundance at the same rate. and the lower atmosphere, there is a rapid accumulation or growth of large-size scatterers. 4. We modeled the phase functions produced by a monodisperse, bimodal, trimodal, log-normal, and power-law distribution of particles. Our results confirmed the empirical analysis of the phase curves, i.e., that a distribution of particles, characterized by an abundance of small size scatterers R f < 50 nm(99%), and a non-negligible amount of large scatterers (100 nm < R F < 1 μm) (1%), which specifically include a crucial particle size of 500-700 nm, are able to fit the phase curves observed below 200 km. Above 200 km, these populations all approximated a monodisperse population, with a characteristic particle size, where Rayleigh scattering dominates type. This can be accomplished by increasing the charge-to-radius ratio to prevent the removal of smaller scatterers, while reducing the sedimentation velocity to allow these more coagulationresistant particles an opportunity to form large scatterers. The sedimentation velocity can be reduced for larger and larger particles with a fractal dimension of less than two.
This work provides constraints from New Horizons data on the population type, abundances, and some physical properties of Pluto's aerosols. A number of avenues for followup are possible, including the following: 1. A full-phase curve from 0°to 180°is needed to determine the abundance of ( ( ) -) » R Int 100 500 nm f , and to break the degeneracy between a power-law, bimodal, and trimodal distribution. This could be accomplished with the inclusion of images from the panchromatic filters or other instruments (i.e., LORRI), or more flyby missions, but an orbiter would provide a more comprehensive and definitive data set. 2. Microphysical models, such as CARMA, will need to consider these phase curves and distributions in future work. Results from rigorous modeling will help us to better understand the processes leading to a bimodal or power-law distribution. Future modeling will also explore the total charge that can build up on these particles, and how they pass their charge on to other aerosol or gaseous components. 3. Tests of the Tomasko et al. (2008) fractal aggregate phase function for a fractal dimension not equal to two have not been verified or reviewed. Validation of other fractal dimensions should be conducted so as to increase confidence in these results.
The raw and calibrated New Horizon data are available via the PDS (Stern et al. 2007). The entire data set can be made available upon request to the corresponding author. The MATLAB implementation of the phase function employed by Tomasko et al. (2008) to determine scattering from fractal aggregates can also be made available upon request to the corresponding author.

Appendix A Data Processing
All of the images are rendered in Flexible Image Transport System (FITS) format, and are readable using standard FITS viewers and software libraries. (Stern 2018). For this work, the raw data is processed into level 2 calibrated products by removing the bias and flat-field pattern to convert pixel values from raw data number (DN) to calibrated DN (Peterson et al. 2017). Level 2 calibrated products were provided directly from coauthor Buratti, and were produced in September of 2017. Radiometric calibration coefficients, stored as keywords in the FITS header, are used to convert raw DN to physical reflectance values. Error estimates for each pixel are calculated assuming shot noise and measured detector parameters (e.g., reading noise and Poisson noise from dark current), and an error array is constructed in a new extension of the FITS file. Finally, the SOC pipeline adds a data quality extension, based on both scene complexity and instrumental factors. The data quality flag is set to zero if the data are clean, and a nonzero value that represents potential issues. These include defects in reference calibration files, such as dead pixels, missing data, zero-value pixels, DN not linear in the detector regime, or an unidentified error (Peterson et al. 2017).
We convert from the calibrated DN to physical reflectance units following the procedures described in C.J. Howett et al. (2017). Navigation files were generated for each image, using the Jet Propulsion Laboratory Navigation and Ancillary Information Facility's SPICE system/routines (Acton et al. 2018). The SPICE system combines kernel files containing information on the spacecraft, camera orientation, planetary ephemerides, etc., with keywords in each image's FITS files, such as time of observation, image size, and scan rate. Using the SPICE API, we determined which pixels should be on the disk, as well as characteristics of the pixels on and off the disk, such as geographic coordinates, tangential altitude, and viewing geometry (i, e, θ, ψ). For pixels off the disk, we also determined the viewing geometry and geographic location of the subtangent point.
Due to pointing and absolute timing uncertainties, the location of our disks, as predicted by SPICE, and the actual location of Pluto's disk seen in the images were shifted by 1-10 pixels in line and/or sample. The same shift was observed for both Charon and Pluto in images observing both simultaneously, suggesting that there is no internal timing error or rotational offset. To account for these random shifts, we manually adjusted the navigation files until the disk predicted by SPICE matched the disk in the image to within 1 pixel radii.

Appendix B Off-disk Glow Correction
There is an instrument effect, similar to the stray light seen in the LORRI (Cheng et al. 2017), which causes spurious brightness far from the limb of the disk of Pluto. This effect can be seen more than 100 pixels away from the disk, equating to more than 500 km in some images. To account for this glow, we follow a similar strategy as that employed by Cheng et al. (2017), whereby we use the airless (i.e., presumably hazeless) body of Charon (Stern et al. 2017) to model the effect. For each observation of Charon in the time frame of our investigation (Table 1), we locate each pixel on the solid limb, and find which pixels off the disk are closest to each of those limb pixels, as shown in Figure 14. Around the disk, we now have brightness profiles, dependent on solid-limb brightness and distance. The solid-limb brightness, I/F L , is taken to be the weighted average of the nearest on-disk pixels to the solid-limb pixel within a radius of four pixels. These weights are determined based on the Gaussian profile of the point-spread function (Reuter et al. 2008). The distances of the pixels from the solid limb, D L , is determined by the altitude of the center of the pixel, divided by the difference in altitude between the highest altitude corner of the pixel, and the lowest altitude corner of the pixel. Both Pluto and Charon cast shadows at photometric longitudes, Λ, for |Λ| > 90°. This means that pixels with a subtangent point of i > 90°, will have a LOS through the atmosphere that is only partially illuminated. To ensure that we are only using pixels whose LOS is totally illuminated, we limit ourselves to those with a subtangent point where i < 90°.
We fit the Charon off-disk brightness profile to be a twodimensional function of the form: Figure 3 highlights our fit to the data, and provides a comparison with Pluto's off-disk profile. Best-fit coefficients, and their associated errors, can be found in Table 2. We find that the magnitude of the removed off-disk glow is on the order of the variability we see in I/F, constituting a substantial fraction of the observed signal. We limit our data for generating phase curves to observations of Pluto where Pluto's weighted limb brightness, I/F L , is within the range of Charon's I/F L . This procedure maximizes the efficacy of our fits, although our functional form of Charon's off-disk profiles allows us to apply this correction to any limb brightness value. We were not able to limit the data in this way for the blue filter, because Pluto is much brighter in the blue than Charon. We also found that variations in the geographic location of the subtangent point coincided with brightness variations at all altitudes, ranging from a factor of two to an order of magnitude ( Figure 15). To limit the variation in the brightness profiles, we constrained our data to subtangent points with an incidence angle of 75°-85°, and in the region of 35°-50°N latitude and 80°W-40°E longitude ( Figure 16).
Phase curves were generated by binning the data into altitude bins 20km high. This ensures that we have at least several pixels in each altitude bin, even at our lowest spatial resolution. We only use pixels that have all four corners within the defined altitude bins. We choose our upper limit to be 500 km, as this is the height where photolysis of methane takes place (Gladstone et al. 2019), and below this altitude, aerosol growth begins (Gladstone & Young 2019). Due to the nature of our Mie scattering code (Mätzler 2002) we must round the phase angle of each pixel to the nearest whole integer. At closest approach, this results in two phase angles, at 40°a nd 41°, so that each of our phase curves have a total of eight points. In the figures presented here, the I/F is presented by the median of the I/F values of the pixels within the altitude bins, and at the same phase angle. For our fitting routines, however, we use a Monte Carlo simulation, with at least 320 runs, to draw our phase curve values from the distribution of pixels in each of the altitude and phase bins.
Negative pixel values arise from level 1 calibrations (i.e., flatfield corrections and dark current subtractions), and the Charon offdisk subtraction. This means that some points in our phase curves have median I/F values below zero, particularly in the NIR channel, which has lowest S/N of our t filters. We excluded these observations from our phase curves, as their inclusion interfered with our fitting routines. Values that were excluded are noted in red in Figures 17 and 18. Recall, however, that this is for the phase curve figures; for our fitting routines we use a 320-run Monte Carlo simulation, in order to draw our phase curve values from the distribution of pixels in each of the altitude and phase bins.   Table 1). (Right) Pixels taken from MVIC sequence 0299193157 (See Table 1    Note. R f : fractal aggregate the effective radius, Equation (15); S f : scaling factor, Equation (7); N: the particle number density Equation (D8); and R 2 : goodness of fit metric, Equation (4). Note. R f , fractal aggregate the effective radius, Equation (15); w: relative weight of aerosol population, Equation (16); S F : scaling factor, Equation (7); N: the particle number density Equation (D8); R 2 : goodness of fit metric, Equation (4).   Note. b: power law exponent; S F : scaling factor, Equation (7); n: the differential particle radius number density distribution, Equation (10); R 2 : goodness of fit metric, Equation (4). Figure 17. Median level 1 calibrated products, converted from calibrated DN to physical reflectance units, I/F, following the keyword procedure detailed in Howett et al. (2017), Binned by altitude and phase angle. Left: data without correction applied from the Charon profile. Right: data with correction applied from the Charon profile applied. Red data is neglected.
In addition to the primary data, each FITS files contains two extensions. The first extension is the 1σ standard deviation of each pixel in DN. The second extension is the data quality flags. These contain the error array, and the data quality array for the primary data. A full explanation of the various flags and error calculation is detailed in New Horizons SOC to Instrument Pipeline ICD. Throughout this work, we only include pixels where the data are flagged as good (Peterson et al. 2017).
Along with the random pixel error, a 1%-3% error exists in the adjustment factor required to correct the observed count rates to the predicted count rates (Howett et al. 2017). This adjustment factor was determined by means of a comparison of photometry with well-calibrated stars observed in-flight by all 7 MVIC CCDs. To calculate error, we combine the Poisson photon noise for each pixel (uncertainty from dark current and bias is negligible) with a systematic multiple offset to account for each frame's absolute calibration error. The calibration error is approximated by a Gaussian distribution, with a standard deviation of 3%.
We present the uncertainty in I/F as the bounds occurring at the 15th and 85th percentile of the pixels in each bin ( Figure 18). The 15th and 85th percentile range is equivalent to the bounds at one standard deviation, σ, greater and less than the mean or median for a normal distribution. The variance in our data in a single bin is typically not normally distributed, so the standard deviation is not appropriate for describing out variance. This variability in our binned data, ∼25% and dependent on filter (Figures 17 and 18), is substantially greater than the 3% propagated error above. As such, we do not consider the propagated error, but we present our determination of the propagated error above for the sake of completeness.

Appendix D Radiative Transfer
A full radiative transfer description of scattering from Pluto's limb would include ray tracing through spherical shells, and account for multiple scattering (Bourassa et al. 2008 which specifies the radiance, I, observed at location, r 0 , in the propagation direction specified by the unit vector,  W, along the LOS path, s. This is the integral of the scattered source term, J, plus the radiance at the end of the path,Ĩ , both attenuated along the path to the observer by the optical depth, τ.Ĩ is total radiance (direct, single scattered, multiple scattered) from the end point of the path, which is either space or the ground (Figure 1 bottom, A and D). As none of our observations intersect the surface of Pluto or the Sun, this term is zero, and will not be discussed further.
The source term from scattering is where β s is the scattering coefficient, P 11 is the scattering phase function (Equation (2)), and θ is the angle formed between incident and emergent vectors, specified by the unit vectors  W¢ and  W, respectively. J is the total scattering source term, where the final scattering event occurs in the atmosphere. The integral is evaluated over all solid angles, as multiple scattering can have its source from any direction, including the Sun (Figure 1(b), B), the ground (Figure 1(b), E), and from all other scatterers in the atmosphere (Figure 1(b), C). Because the solar beam is the only source of light for the first-order scattering, the integral over all space that is required to evaluate the source term (Equation (5) where ( )   W F 0 is the incident solar flux density. Thus the total contribution to the radiance at any position r 0 from a single scattering event is Pluto's, however, the integral in Equation (4) is largely dominated by scattering from the tangent point, reducing Equation (4) to Normand et al. (2013): where s T , is the tangent point in Pluto's atmosphere (for our purposes, s T is the altitude), and Δs T is the tangent point path length. Accounting for all components of Pluto's atmosphere Equation (10)  where λ is the pivot wavelength of our MVIC filters, N(s T ) is the particle number density at s T in units of per volume, and σ is the scattering cross-section in units of area. The particle number density for Pluto's atmosphere is described in Young et al. (2018), based on a UV solar occultation from the New Horizons Alice instrument to measure the abundance of N 2 , CH 4 , C 2 H 2 , C 2 H 4 , C 2 H 6 , and haze, as a function of altitude. They show that abundance in Pluto's atmosphere is dominated by N 2 , with 10 15 cm −3 , and CH 4 , with 10 12 cm −3 at the surface. We determined the magnitude of the scattering due to Pluto's gaseous components using these abundance altitude profiles, a scattering crosssection for N 2 and CH 4 from experimental work ( We determined the magnitude at 370 nm, where laboratory data for N 2 and CH 4 cross-sections exists. This is 122 nm shorter than our BLUE filter pivot wavelength, and according to the Rayleigh scattering intensity approximation of I ∝ 1/λ 4 , the intensity at 370 nm is 3 times stronger than at 492 nm. Even at this wavelength, we find that the Rayleigh scattering at the surface (0 km) due to N 2 is I/F ≈ 10 −3 , while Rayleigh scattering due to CH 4 is I/F ≈ 10 −7 . The mean variability for the MVIC Blue filter at 0-20 km for all phase angles is ≈10 −2 Given that these values for the UV wavelength, and at the surface where n is at its largest, are at most an order of magnitude lower than the variability we see with the MVIC Blue Filter, we can ignore the contribution from the gaseous components. This is even more the case for the Red and NIR filters, given that I ∝ 1/λ 4 , i.e., the contribution from the gaseous components, is almost nonexistent. Therefore, Equation (11)

Appendix E Phase Function
Many analytic and numerical techniques have been developed to solve for the scattering of electromagnetic waves by arbitrary scatterers (Sharma 2018). These exact solutions, however, can be complicated to implement, and may have dubious analytic value when trying to diagnose a scatterer from a phase curve, especially if the phase curve is discontinuous or incomplete.
We investigated many relevant phase functions, including Mie scattering, the Henyey-Greenstein phase function, various modification of both Mie scattering and Henyey-Greenstein that more accurately approximate real scatterers, and models that handle extreme or niche situations (Henyey & Greenstein 1941;Cornette & Shanks 1992;Liu 1994;Draine 2003;Zhao et al. 2006;BenZvi et al. 2007;Aartsen et al. 2013). Following this exploration, we decided to use the empirical phase function described by Tomasko et al. (2008) for scattering by a fractal aggregate composed of small monomers in Titan's atmosphere. Fractal aggregation is consistent with the color images of Pluto, which exhibit strong forward scattering by a predominantly blue atmosphere (Gladstone et al. 2016). Titan's haze seems to be composed of fractal aggregates, formed as a consequence of photolysis, and coagulating as they travel through Titan's atmosphere (e.g., Lavvas et al. 2013;Gladstone et al. 2016;Cheng et al. 2017;Young et al. 2018); given the similarities in atmospheric composition, this suggests that Pluto's haze may share a similar composition.
The Tomasko et al. (2008) phase function uses the fractal dimension, D f , the number of monomers per aggregate, N, the complex index of refraction as input parameters, m r + im i , and the monomer size parameter, X m , given by X m = 2πr m /λ, where r m is the spherical radius of the monomers. The tested range for these parameters can be found in Table A1 of Tomasko et al. (2008). As a part of this work, we constructed a Matlab implementation of the phase function described in the appendix of Tomasko et al. (2008).
The composition, given by proxy through the complex index of refraction, does not change the phase angle brightness dependence, or the shape of our modeled phase curves. It does, however, change the scattering and extinction cross-sections, σ s and σ e , the scattering and extinction coefficients, β s and β e , and the particle number density for particles of a radius r, N(r) (Equation (13)), which all dictate the absolute brightness of our phase curves. As such, we follow the lead of other authors Young et al. 2018) and assume that the haze has the same composition as Titan's haze. We use the complex index of refraction for Titan's haze from Ramirez et al. (2002). The wavelength range covered by Ramirez et al. (2002) ends at 200 nm. To use the abundance measurements made by the Alice instruments UV solar occultations (Section 3.3. Size distribution), it was necessary to linearly extrapolate the complex index of refraction from Ramirez et al. (2002) to 165 nm. We found the real and imaginary components of the Titan aerosol at 165 nm to be m r = 1.6839 and m i = 0.0166.