Spectral Modeling Using Radiative Transfer Theory with Packing Density Correction: Demonstration for Saturnian Icy Satellites

We demonstrate the capabilities of the radiative transfer theory with packed media correction (RTT-PM) in analyzing spectral data of planetary surfaces by modeling to first order the shape and band depths of spectra of icy satellites of Saturn acquired by Cassini Visual and Infrared Mapping Spectrometer (VIMS). The RTT-PM is an efficient and physically strict numerical method that employs a packing density correction, the static structure factor, to single-scattering properties of particles to simulate the light scattering by densely packed media. Originally created for layers formed by spherical homogeneous particles, the RTT-PM method has been recently updated to treat particles of arbitrary shapes and structures, including aggregates. We apply the RTT-PM method to roughly model Cassini VIMS spectra from Dione, Rhea, and Tethys as layers of spherical particles versus aggregates. The shape and structure of particles strongly affect the modeled spectra; the best model comparisons to the VIMS spectra were obtained when the surface icy particles were assumed to be small aggregates consisting of micron-sized monomers, which may imply rather compact, irregular particles. Our results suggest that presenting the icy regolith as a dense layer of nonspherical particles may noticeably affect the modeling results and bring a better understanding of the satellite surface structure and composition. The RTT-PM demonstrated itself to be a powerful tool for such studies: we computed a reflectance for 22 wavelengths within minutes using a regular desktop computer. The combination of such high efficiency and physical strictness makes the RTT-PM method advantageous for analyzing large spaceborne instrument data sets.


Introduction
Interpretation of the light scattered by particulate (regolithic) surfaces is one of the most important and complicated problems in remote sensing. For example, one of the major means to study solar system bodies is by analyzing reflectance or the similar radiance/irradiance ratio (I/F) spectra from spacecraft or Earth-based remote sensing measurements. By modeling reflectance (hereafter includes I/F) spectra, one can characterize compositional and physical properties of planetary surfaces, which is essential in solving various mysteries existing in the solar system and beyond. However, modeling reflectance spectra of particulate surfaces such as regolith is complicated by multiple scattering, a process in which light interacts with many particles before it is measured. Light detected from multiple-scattering media is affected by physical properties of the media, e.g., the packing density and particle size, shape, and structure, in addition to the geochemical composition of the materials.
Multiple scattering is often modeled with a radiative transfer approach (e.g., Chandrasekhar 1960). Classical radiative transfer theory (RTT) describes the multiple scattering of light in a particulate medium with a far-field assumption where individual particles are considered to be well separated enough from each other, as they are, for example, in clouds and colloids (see Mishchenko et al. 2006). The far-field approximation requires that particles be separated by at least 3 particle radii (van de Hulst 1957). This is not the case for planetary surfaces where many particles exist in close proximity or are in contact with each other. Nevertheless, various assumptions and modifications can be made to the traditional RTT to force it to be applicable to planetary surfaces, and radiative transfer modeling of spectra remains an important method in characterization of planetary surfaces.
The radiative transfer models most familiar to planetary researchers are the analytic isotropic multiple-scattering approximation and the anisotropic multiple-scattering approximation models developed and refined by B. Hapke and collaborators (see overview of this approach in Hapke 2012 and references therein). The Hapke approaches allow efficient computation of photometric and spectral characteristics of light scattered by particulate surfaces. However, they involve approximations and parameters that are not the physical properties of the surfaces, which are compositions, sizes, shapes, and structures of particles and of the surface, but rather their derivatives (single-scattering albedo, asymmetry factor, and width and strength of the opposition surge). The derivative properties are not related to the physical surface properties in a linear and unique fashion. For example, the single-scattering albedo and packing density are mathematically coupled and can produce nonunique solutions (Shepard & Helfenstein 2007). This complicates or even makes the extraction of the desired physical surface properties from Hapke modeling impossible.
Another popular approach that is not based on RTT is the superposition T-matrix approach (MSTM; Mackowski & Mishchenko 1996Mackowski 2014). It is based on a solution of Maxwell's equations and, thus, rigorously includes all phenomena involved in light scattering by media. It is directly based on physical characteristics of the media: the particle size, the particle complex refractive index, and the media's packing density. However, the superposition T-matrix approach is very demanding on computational resources when modeling large volumes of particulate media or media consisting of large particles. Even in its parallelized version, the MSTM code requires a large amount of computer memory and wall clock time for computations to complete and thus can be practically unfeasible for analyzing large data sets. Its current version is also only valid for spherically shaped (although polydisperse and inhomogeneous) particles.
Note that there is also the radiative transfer coherent backscattering (RT-CB) technique (Muinonen 2004;Muinonen & Videen 2012;, which is a Monte Carlo radiative transfer modeling solution with a control of the phase difference between the waves propagating along the reciprocal paths in the medium. It has been used so far mainly to interpret the light scattering at small phase angles (opposition phenomena). Besides, it works only for sparse media, and it is also demanding on computer resources. One more limitation of the RT-CB approach is that, unlike in the Hapke or MSTM cases, its codes have not been publicly available until very recently, constraining the applications of the approach by its developers.
A great compromise between the efficiency and rigor of the approaches considered above for analyzing spectral measurements is the approach presented in this work: radiative transfer theory with packed media correction (RTT-PM). The RTT-PM combines a radiative transfer solution with a rigorous packing density correction and is applicable to both spherical and more complex shaped particles. The basic codes used in this technique are available online at https://www.giss.nasa.gov/ staff/mmishchenko/brf/. This approach has been successfully used to study light scattering by densely packed media such as mineral powders and snow (Tsang et al. 1985;Tsang 1992;Mishchenko 1994;Pitman et al. 2005;Ito et al. 2018). In this paper, mainly directed to show the capabilities of the method, we use the RTT-PM method to create first-order models of the Cassini Visual and Infrared Mapping Spectrometer (VIMS) spectra of three high albedo icy satellites of Saturn. The main purpose of this exercise is to see how shape and structure of icy particles affects the spectra; specifically, we are trying to find out if the icy regolith of those objects is formed by solid icy grains or by aggregates, which are typical for comets and some other icy bodies.
In Section 2, we describe the RTT-PM method and outline the modeling procedure. In Section 3, we describe the observational data and compare them with the results from RTT-PM modeling in Section 4, where we also discuss our findings. Section 5 summarizes the conclusions of the paper.

Brief Theoretical Background
The RTT-PM method is based on the phenomenological RTT (Chandrasekhar 1960). In short, the phenomenological RTT accomplishes the derivation of the scattering properties of the whole particulate medium by integrating the scattering characteristics of fundamental particles constituting the medium, which are namely the single-scattering albedo (the ratio of scattering and extinction cross sections) and the 4×4 normalized Stokes 5 scattering matrix. A number of wellestablished methods are currently known to compute the singlescattering albedo and the scattering matrix for single particles such as, but not limited to, Mie and the superposition T-matrix methods used in this work. Numerous other techniques to compute these quantities, e.g., discrete dipole approximation (DDA; Draine & Flatau 1994, see updated code versions athttp://ddscat.wikidot.com/), are described in the review by Wriedt (2009). The phenomenological RTT is traditionally a theory for sparsely packed particulate media, meaning that light scattered by one particle is assumed to not interfere with other particles existing in the medium (independent scattering assumption). The single-scattering albedo and scattering matrix of a particle in the RTT are given in the far-field zone by definition, but these far-field quantities are not strictly appropriate for particles in a densely packed medium where they exist in close proximity to one another. Notably, the forward scattering intensity of a particle in the independent scattering assumption is too large whereas dense packing leads to lower scattering in the forward direction (e.g., Mishchenko 1994), and thus, the use of the single-scattering properties derived from common methods like Mie and the T-matrix techniques are not entirely appropriate for RTT in densely packed media.
To overcome this problem, the RTT-PM method modifies the single-scattering parameters to include the interdependence effects of closely located particles using a technique known as the static structure factor correction (see mathematical derivations in Tsang et al. 1985;Mishchenko 1994;Pitman et al. 2005;Ito et al. 2018). The interaction of particles with one another is modeled using the Percus-Yevick approximation of the pair-distribution function (Tsang et al. 1985;Tsang 1992;Mishchenko 1994). This is a statistical mechanics-based approach where a radial distribution function (Percus-Yevick pair-distribution function) is used to get the probability of finding a particle at a given distance from a reference particle. For a certain wavelength and particle size, the probability of finding a particle (thus the structure factor) as a function of scattering angle depends on the packing density of the medium, i.e., the volume fraction of the medium occupied by the particles (often called the "filling factor" in the previous papers). Physically, the structure factor accounts for the interference of light scattered by the particles such that more densely packed media lead to stronger interference effects, and it redistributes the scattered energy, which leads, for example, to a reduction in the forward scattering peak.

Computational Technique
The RTT-PM method computes reflectance in a three-step procedure: 1. Calculate the single-scattering parameters, i.e., the singlescattering albedo (from scattering and extinction cross sections) and Mueller scattering matrix elements, of a single particle that can be solid or represents a cluster of grains (aggregate).
2. Apply the static structure factor correction to the singlescattering quantities from step 1.
3. Compute the bidirectional reflection function (or hemispherical averages) using a radiative transfer model with the results of steps (1) and (2).
In this study, for step (1), we use both Mie and the superposition T-matrix methods. The computer codes for these methods are available online at https://www.giss.nasa.gov/ staff/mmishchenko/ftpcode/spher.fandhttp://www.eng.auburn. edu/∼dmckwski/scatcodes/, respectively. Mishchenko (1994) used the RTT-PM approach with Mie modeling as the first step to compute scattering by snow and soil. Mie theory (Mie 1908;Bohren & Huffman 1983) describes the characteristics of light scattered from a homogeneous spherical particle using only physical characteristics (the wavelength of light, complex index of refraction for the given wavelength, and particle radius) as model inputs. Ito et al. (2018) extended the step (1) approach to nonspherical particles, specifically to clusters of closely spaced spheres (aggregates) whose scattering properties can be rigorously calculated using T-matrix methods by Mishchenko et al. (2002) and Mackowski & Mishchenko (1996. The superposition T-matrix method known as Multiple Sphere T-Matrix (MSTM) rigorously solves Maxwell's equations for the cluster of particles as a whole, leading to incorporation of all effects, including near-field effects, which can be crucial to have when modeling light scattering by densely packed media (e.g., Tishkovets 2008; Ito et al. 2017). The inputs to the T-matrix codes are again the wavelength of light, complex refractive index for the given wavelength, radius of the constituent sphere, number of spheres in the cluster, and packing density of the cluster (i.e., the ratio of the volume of constituent spheres to the total volume of the cluster), usually defined through the porosity (i.e., unity minus packing density). In this paper, we use the terms "packing density" for the medium and "porosity" for aggregates. The outputs from these codes are the single-scattering albedo (as well as the scattering, extinction, and absorption cross sections necessary to derive it), the 4 × 4 Mueller scattering matrix, and a transformed version of the scattering matrix known as generalized spherical function (GSF) expansion coefficients.
To model realistic surfaces of planetary moons and asteroids, we need to consider media consisting of rather densely packed particles. In step (2), we apply the static structure factor correction, described in Section 2.1, to the single-scattering albedo and scattering matrix obtained in step (1). The second step corrects the phase function. This causes a reduction in the large forwardscattering peak that automatically appears as part of the Mie or T-matrix methods due to these techniques reporting output scattering quantities in the far-field zone by definition.
The static structure factor correction is accomplished with the computer code from Ito et al. (2018; available on https://www. giss.nasa.gov/staff/mmishchenko/ftpcode/PackedMedia.f). The inputs to this code are the cross sections (scattering, extinction, absorption), GSF expansion coefficients, and the number of GSF expansion coefficients, which are all provided as the results of step 1. The packing density should also be specified. The outputs are the single-scattering albedo, scattering matrix, and GSF expansion coefficients that have been corrected by the static structure factor.
The third step of the RTT-PM method is the computation of reflectance (the bidirectional reflection function) using the radiative transfer equation. Using the single-scattering properties of the individual particles (step 1) with dense packing modifications (step 2), the radiative transfer (step 3) calculates light scattering by the entire medium. In this study, we assume that the surface under consideration is macroscopically flat, homogeneous, and optically semi-infinite, and we apply the invariant-imbedding solution to the scalar radiative transfer equation (Mishchenko et al. 1999). The invariant-imbedding RTT method considers the change in reflection from a medium by the addition of an optically thin layer in which the thin layer scatters light at most once (summarized in Hansen & Travis 1974;Mishchenko 2014). The solution is obtained numerically by iterating the Fourier-decomposed Ambartsumian nonlinear integral equation as shown in Mishchenko et al. (1999). The computer code for this step is available on https://www.giss.nasa.gov/staff/mmishchenko/brf/ where refl.f outputs the Fourier components of the reflection function and interp.f transforms the output from refl.f into the bidirectional reflection function. Hemispherically averaged reflectance is also included in the output from refl.f, although we do not use it in this study. A detailed manual can be found in Mishchenko et al. (1999Mishchenko et al. ( , 2015. Figure 1 outlines the master plan of the RTT-PM computations with the specific parameters used in this paper noted. Thus, the RTT-PM procedure uses the complex refractive index, particle radius (or the radius of constituent monomers and their number and porosity for aggregates if using the superposition T-matrix method), and packing density of the medium to model reflectance spectra of planetary surfaces. The usage of a few key physically based parameters important for modeling reflectance spectra of regolithic surfaces and the method's reasonable computational practicality are primary advantages of the RTT-PM method.

Surfaces of Midsize Saturnian Icy Satellites Studied by the Cassini VIMS Instrument
In this study we apply the RTT-PM technique to simulate the shape and band depths of the near-infrared (NIR) spectra of icy satellites of Saturn. These objects are especially appealing for our study because the RTT-PM method has not yet been validated for outer solar system data. Also, we are testing RTT-PM on a few VIMS spectra previously studied with the MSTM modeling technique in the hope that we can produce similar results (i.e., preserve the physics but at a lesser computational investment) and therefore contribute an additional tool for the outer planet community to analyze large spacecraft data sets.
We apply the RTT-PM method to create model reflectance spectra of the midsized icy satellites of Saturn: Rhea, Dione, and Tethys, which are rich objects of study that offer opportunities to compare and contrast different surface and structural properties to better understand geology, formation history, and evolution of moons. Buratti et al. (2018) presented an excellent summary of the open problems to be solved for the Saturnian icy satellites. Rhea, Dione, and Tethys are often grouped together in studies because they all orbit inside of Saturn's E-ring, giving rise to leading hemispheres that are brighter in the visible wavelength range (McCord et al. 1971;Buratti et al. 1990;Verbiscer & Veverka 1992). They also have similar systems of faults and craters (Beddingfield et al. 2015(Beddingfield et al. , 2016Byrne et al. 2015;Neidhart et al. 2017), and their particular masses and orbital inclinations have implications in models of satellite formation (e.g., Ćuk et al. 2016;Salmon & Canup 2017).
We use Rhea, Dione, and Tethys as convenient samples to test the RTT-PM technique on real-world planetary spectra. There are two main reasons for this. First, their spectra were the subject of modeling in previous papers (Kolokolova et al. 2010(Kolokolova et al. , 2011Pitman et al. 2017), where the MSTM techniques, specifically the MSTM3 code for large volumes filled with spheres (Mackowski & Mishchenko 2011) and the MSTM4 code for layers of spheres (Mackowski 2014), were used. Thus, with this study we can compare the results obtained by the rigorous MSTM solution and RTT-PM technique to see if a much more computationally efficient RTT-PM approach accounts for the major light-scattering phenomena that are responsible for the formation of the reflection spectra of the icy satellite surfaces. Second, those three satellites do not have a significant atmosphere that could complicate the modeling. Although Rhea and Dione were found to have tenuous seasonal exospheres composed of oxygen and carbon dioxide concentrated over their dayside hemispheres, as confirmed by Cassini instruments Cassini Plasma Spectrometer (CAPS) and Ion and Neutral Mass Spectrometer (INMS; Teolis & Waite 2016 and references therein), their contribution should be negligible for this application, especially if we consider the high spatial resolution VIMS data taken from close distances and covering small surface areas. In this paper, we focus on a single composition (water ice) while varying structural characteristics of the particles and the medium to compare against previous work and to make the modeling more straightforward. The relatively high albedo of the three icy satellites and prominent bands of water ice in their NIR spectra allow us to suggest that they are composed of predominantly water ice. The full composition of icy satellites is believed to include contaminant species such as carbon, tholin, complex organics, amorphous silicates, nanophase iron, and hematite (see Johnson et al. 1983;Moore et al. 1983;Cuzzi & Estrada 1998;Poulet et al. 2003;Clark et al. 2008;Ciarniello et al. 2011), which may affect albedo and the general slope of the spectra, as well as cause "dents" on the slopes of some bands for those modeling compositions in detail. Modeling results by other groups (e.g., Clark et al. 2008;Ciarniello et al. 2011;Filacchione et al. 2012) suggest that the amount of these contaminants is small, at the 0.4%-0.5% level, which does not have a large impact on our analysis here. The Dione and Tethys VIMS spectra shown in Pitman et al. (2017) were acquired on the trailing side of those moons. Dione's trailing hemisphere ranges from contaminated to pure water ice compositions. Of the spectra available from the Pitman et al. (2017) study, we have selected a Dione spectrum extracted near a pure water ice region (VIMS spectral unit # 8; Scipioni et al. 2013) to minimize the impact of contaminants. The degree of contaminants in the Tethys spectrum from that study is unknown. The Rhea spectrum from that study has some degree of contaminant (likely VIMS spectral unit # 2 or 3 composition; Scipioni et al. 2014). The effects of contamination can be easily included in the RTT-PM modeling by modifying the refractive index or by considering a mixture of different particles and is planned to be considered in future work. However, in this paper, we avoided adding contaminants to achieve a spectral fit as this could mask the potential contributions from particle shape and structure. Besides, focusing on the icy surface allows us to compare our RTT-PM modeling with the previous work based on the rigorous MSTM approach (Kolokolova et al. 2011;Pitman et al. 2017).  Table 1 lists the Cassini VIMS observations of Rhea, Dione, and Tethys used as the "bench spectra" to test RTT-PM. Details of the data selection and processing were described in Pitman et al. (2017). Briefly, the Cassini VIMS instrument acquired a large hyperspectral data set over a spectral range of λ=0.34-5.1 μm (Brown et al. 2004), with a 256 band IR channel (0.88-5.1 μm), a spectral resolution of Δλ ≈ 16 nm/ band, and spatial resolutions of 0.5 × 0.5 (nominal) or 0.25×0.5 (high-resolution) mrad/pixel (Miller et al. 1996). The spectra for our modeling were extracted from Cassini VIMS data cubes archived at the NASA Planetary Data System (PDS 2016), being careful to avoid cubes with missing pixels and low signal-to-noise ratios. The VIMS data were calibrated using U.S. Geological Survey Integrated Software for Imagers and Spectrometers (ISIS 3), translating the offset between the geometric backplane information and the icy satellite image as instructed by the Cassini VIMS science team. The specific spectra from Pitman et al. (2017) selected for this study were chosen to allow a direct comparison between RTT-PM and MSTM4 numerical techniques. The main advantage of those VIMS spectra is that they are not from the whole area bounded by the VIMS observation but rather for specific pixels within it (Pitman et al. 2017, see red or blue dots in Figure 3 there); thus, they were acquired for rather small areas on the icy satellites' surfaces that can be considered as flat. Also, the icy satellite surfaces in those areas do not cover large-scale geologic features, making them suitable for modeling a plane-parallel medium. The geometrical characteristics of the considered surfaces and of their observations are listed in Table 1.
We selected the observations at solar phase angles close to 1°to allow for further comparison with the results of Pitman et al. (2017) where this was the only similar solar phase angle for which we had the data for all three satellites. Besides, for icy satellites, the opposition surge caused by the coherent backscattering occurs at a solar phase angle <1° (Buratti et al. 2009). In order to minimize the variability induced by the solar phase, we have selected VIMS observations taken at very similar angles (Table 1).

Model Inputs
The RTT-PM technique described in Section 2 was used to model the Table 1 VIMS spectra. The main purpose of our modeling was to use the unique capability of the RTT-PM to model media formed by particles of different (but controlled) shapes and structures. In order to focus on the effects of shape and structure of particles, we limited the modeling to monodisperse particles made of pure ice. We started our modeling by assuming the icy satellite surfaces were a medium formed by spherical particles made of water ice. Following the benchmark modeling in Pitman et al. (2017), the wavelength-dependent refractive indices of water ice were taken from Mastrapa et al. (2008Mastrapa et al. ( , 2009. We simulated the spectra for spherical particles with radii from 0.5 to 5 μm using a 0.5 μm step in radius and found that the match to the VIMS spectra improved when we did the calculations for the sphere with a radius of ∼2 μm. Then we repeated this procedure with a smaller step in the radius to find the best representative size. However, we know that the regolith particles are usually nonspherical and may have a complicated shape and structure. In many cases, they are aggregated (agglomerated) particles; we know this, for example, from the studies of interplanetary dust particles of asteroidal and cometary origin (e.g., Sandford & Bradley 1989;Güttler et al. 2019) and laboratory simulations of regolith light-scattering characteristics (e.g., Piatek et al. 2004). Thus, we also investigated whether the surface of the icy satellites is covered by aggregated particles.
We represented the aggregated particles as clusters of spheres and built them using a ballistic approach, considering ballistic particle-cluster aggregate (BPCA) and ballistic cluster-cluster aggregate (BCCA) models (Meakin 1984), resulting in aggregates of porosity of about 85% and 95%. The light-scattering properties of a single aggregate were computed using the T-matrix MSTM code for randomly oriented clusters of spheres, which has been made available online at http://www.eng.auburn.edu/ dmckwski/scatcodes/ by its developer D. Mackowski. Similar to the procedure taken for the single-sphere case, we considered aggregates formed by monomers (constituent spheres) of radii from 0.5 to 5 μm with the step in the constituent sphere radius of 0.5 μm. We started with large aggregates (1024 spheres in a cluster) but found that the results were very different from the observed spectra; absorption bands were either too shallow or even absent. Then we considered the smallest clusters of two spheres and found that with them we can achieve a much better match to the measured spectra. After that, we started increasing the number of constituent spheres in the aggregates and found that the model comparison to the VIMS observed spectra usually became worse as the number of spheres in the aggregates increased, except for the case of Dione where aggregates of 128 constituents led to the best model to real-world data comparison. Larger aggregates did not work for all satellites.  parameters. The quality of the fit is presented by root mean square (rms) error calculated as a square root of the sum of squared differences between measured and calculated values of normalized reflectance for each of 22 wavelengths, which form the spectra shown in Figure 3, divided by the number of wavelengths.

Results of the RTT-PM Modeling
The results for aggregates with the rms error exceeding unity were not included in the plot; these were usually the aggregates of a large monomer radius (>2 μm) and/or a larger number of monomers in the aggregates (usually >128). Note that in the case of spheres, a smaller packing density of the medium produces better results; the same trend was found for aggregates. It can be seen from Figure 2 that the change in the rms error for aggregates was not smooth when the parameters of the model smoothly increased or smoothly decreased. This may indicate that not only size of monomers and their number but also the aggregate shape (that is different in the case of different size/number of monomers) affects the results. The best model results for both spheres and aggregates are presented in Table 2 and shown in Figure 3. Figure 3 provides some additional information about the fit presenting the residuals between measured and modeled data for different packing densities.
Note that for the MSTM4 runs in Pitman et al. (2017), a single MSTM4 run took 50 hr and the waiting time in the NASA Advanced Supercomputing (NAS) facility queue was on average six days (sometimes up to two weeks), so the focus for that study was on a single, 2 μm absorption band. The RTT-PM computations in this study were performed for the 1.2-3.2 μm range of wavelengths that covered three absorption bands. Using a simple Dell desktop with Intel(R) W5590 @ 3.33 GHz CPU and 4 GB RAM, all steps of computations together took no more than 10 minutes for a single spectrum. Therefore, the RTT-PM method could be applied to a wider range of scenarios with varying compositions while maintaining a quality comparable to the rigorous MSTM approach and expanding our capability of analyzing Cassini VIMS or other remote sensing data.

Discussion
Our best model results are close to those obtained with rigorous MSTM4 solution for the same satellites of Saturn (Pitman et al. 2017) where a surface was modeled as a layer of spherical particles and the radius of particles was found to be    Table 2 for model input parameters). Different colors indicate different packing densities of the medium.
2 μm with porosity of 95% for Rhea and Tethys and 85% for Dione. Thus, the computations performed with the less computationally intensive RTT-PM method resulted in similar results to those produced by a rigorous solution of the Maxwell equations where all effects of electromagnetic interaction between the particles in the media were taken into account. This is good evidence of the physical robustness of the RTT-PM technique. The advantages of the RTT-PM approach is not only its high computational efficiency discussed above but also its capability to consider particles of complex shape and structure that is not possible with the MSTM technique, which requires the regolith layer to be formed by homogeneous spherical particles. As described in Section 4.1, the main goal of this study was to check the capability of the method to model media formed by particles of different shapes and structures and to verify that it could reproduce the general spectral properties of planetary surfaces, specifically for the outer solar system. For this, we considered two extreme cases: the medium formed by spherical homogeneous particles and by aggregates of different sizes and the number of monomers and different porosity. To stress the effect of particle shape, we used a simple model of the regolith on icy satellites of Saturn, presenting it as a medium formed by monodisperse spherical or aggregated particles made of pure water ice. The main result of the modeling is that the particles on the surface of Dione, Rhea, and Tethys are not large aggregates as would be typical, for example, for cometary icy particles (e.g., Kelley et al. 2013). The models that best approximated the real-world observations appeared to be small aggregates, consisting of 2-4 spheres.
An interesting observation is that with MSTM4, Pitman et al. (2017) could not obtain a modeled spectrum for Dione with rms error values as low as those for Rhea and Tethys when modeling with spherical particles. In this study, we also found that the model for Dione is very different from the model for Rhea and Tethys. This confirms that Dione's surface has some noticeable light-scattering differences from the surfaces of Rhea and Tethys in the areas studied. Pitman et al. (2017) explained specifics of the modeling for Dione by different macrostructure of Dione's surface that could be seen in the image of the studied area. The spectral extraction for the Dione VIMS data was done at a site adjacent to a pure ice region to minimize the influence of contaminants; however, there may still be some contribution from a non-ice material. The exact spot occurs under a text label on Scipioni et al. (2013)ʼs Figure 8; therefore it is unclear whether this specific site corresponds to the neighboring pure ice VIMS spectral unit (#8 dark green in Scipioni et al. 2013ʼs classification scheme) or contains a degree of contaminant. The scattering on Dione requires a larger absorption of light, which may result either from presence of a dark material or from multiple scatterings in larger aggregates or by the surface roughness.
Our results allow us to characterize the particles of icy regolith as compact irregular particles. This is important not only as new knowledge about the icy regolith particles but also because the computations for small clusters of spheres or for simple nonspherical particles are CPU time and memory efficient and can be easily performed using nonparallelized versions of DDA, MSTM3 or other T-matrix codes and regular desktop computers. Combining computer efficiency of the RTT-PM method with the low-demand single-particle computations, the model can be extended by considering polydisperse particles and/or a medium of more complex composition that includes such components as carbon, organics iron, silicates, and other contaminants suggested in Clark et al. (2008), Ciarniello et al. (2011), Filacchione et al. (2012, and other studies. Moreover, the RTT-PM modeling can be done for a mixture of particles of different shape (see Vázquez-Martín et al. 2020) and composition as well as for aggregates made of monomers of different size and composition, or of monomers formed by an intimate mixture of materials. The model can be further advanced by considering particles made of birefringent or optically active materials, which are available as options in the MSTM and DDA codes.

Conclusions
The RTT-PM technique is an efficient and physically strict technique to model light scattering from densely packed particulate surfaces. The input parameters of the model are the size of the particles, their composition, and the packing density of the medium.
The improved method allows modeling of light scattering by media that consist of particles of arbitrary shapes, including aggregates. After the light-scattering properties of a single randomly oriented aggregate are calculated using T-matrix or other techniques, the output scattering parameters of those aggregates is used as input parameters to the static structure factor procedure, which is then used in the RTT. The final output of the current version of the RTT-PM method is the scalar bidirectional reflection function.
Thus, the modern version of the RTT-PM technique has three features that make it particularly well suited to analyze the photometric and spectral data related to densely packed media, especially in the case of large data sets: (1) Efficiency of the technique: the computations can be done using ordinary desktops/laptops within several minutes.
(2) All input parameters represent basic, physically meaningful characteristics of the medium (the composition i.e., optical constants, size distribution, and packing density).
(3) Capability to model media formed by particles of a variety of shapes and structures.
As a test case, we modeled spectra of icy satellites of Saturn using the NIR spectra acquired by Cassini VIMS for specific areas on the surface of Dione, Rhea, and Tethys. When the modeling was done using media formed by spherical particles, we could achieve reasonable approximations to the observed VIMS spectra for icy particles of a radius of 1-2 μm and a small/medium packing density (below 0.4). When we replaced the spherical particles by clusters of spheres (aggregates), only the modeled spectrum of Dione noticeably improved. The best model for Dione was for a 85% porosity aggregate consisting of 128 spheres, each with a 1 μm radius. For Rhea, the best aggregate model was just a bi-sphere with 1.5 μm radius constituents, although its rms error was not so different from the case of a single 2 μm radius sphere. For Tethys, the best model was a four sphere cluster where each sphere had a radius of 1.5 μm. However, unlike the pattern observed for Rhea and Dione, the best rms error provided by the aggregate was higher than the one provided by a case of a single 2 μm radius sphere.
We suspect that the surfaces of Rhea and Tethys are most likely covered by solid icy particles, which more likely have a shape that slightly (in the case of Tethys) or strongly (in the case of Rhea) deviates from the spherical one rather than small aggregates. The RTT-PM models for Dione show that its particles may be small aggregates. However, the Dione spectra can be affected by the surface macrostructure that increases the multiple scattering mimicking the effects produced by multiple scattering in the particles of a complex structure or by the presence of some dark material that causes absorption comparable with that caused by multiple scattering.
In general, our computations showed that the packed layers formed by nonspherical, particularly aggregated, particles are characterized by the spectra that can be significantly different from the spectra of similar layers but formed by spherical particles.
We are planning to continue our modeling considering more complex composition and size distribution of the particles, although our hope is that this study stimulates other researchers to continue similar studies using the RTT-PM for modeling more complicated cases.