Near-infrared Photometry of the Moon’s Surface with Passive Radiometry from the Lunar Orbiter Laser Altimeter (LOLA)

Examining the reflectance of the Moon’s surface across a broad range of viewing geometries through photometric analysis can reveal physical and geological properties of its regolith. Since 2013 December, the Lunar Orbiter Laser Altimeter (LOLA) on board the Lunar Reconnaissance Orbiter (LRO) has been operating as a near-infrared (1064 nm) passive radiometer when its laser is turned off. We present a new analysis of this data set spanning roughly 8 yr and covering the surface up to high latitudes in both hemispheres. We apply semiempirical phase functions to find a lower photometric slope and a narrower opposition effect for the highlands than the maria, consistent with theoretical expectations given the higher albedo of the highlands. Examining various geological properties at global scales shows that, in the highlands, iron abundance (FeO) and optical maturity (OMAT) are the dominant factors affecting the phase function, with a smaller influence from surface slope. In the maria, FeO is the dominant factor, with smaller influences from OMAT, surface slope, and TiO2. Submicroscopic iron abundance (SMFe) has a similar effect to OMAT in both highlands and maria. Analysis at specific sites, including the Reiner Gamma swirl and several silicic anomalies, indicates that the phase functions are consistent with the global data for similar FeO and OMAT. Thermophysical properties inferred from surface temperature observations by the Diviner Lunar Radiometer Experiment on board LRO do not affect the 1064 nm phase function, possibly due to a difference between their depth scale and LOLA’s sensing depth.


Introduction
The reflectance of the lunar surface is a function of wavelength, viewing and illumination geometry, the scattering behavior of the regolith particles, and the microtexture of the surface (Shkuratov et al. 2011;Hapke 2012a).These properties hold vital clues to many surface processes operating on airless bodies throughout the solar system.Therefore, it is crucial to fully characterize the surface reflectance over as wide a range in wavelengths and geometries as possible.To that end, the Lunar Orbiter Laser Altimeter (LOLA) on board the Lunar Reconnaissance Orbiter (LRO) has been making photometric measurements of the lunar surface for over 8 yr.In addition to measuring lunar topography, LOLA can act as a narrowband near-infrared passive radiometer when its laser is turned off.LOLA is now the longest-running such instrument to orbit the Moon, and, as such, it provides a wealth of reflectance data extending LRO's wavelength coverage for photometric studies.Barker et al. (2016) analyzed the first year of measurements, which were restricted to the northern hemisphere.They concluded that iron abundance and optical maturity (OMAT) were the dominant geological factors controlling the nearinfrared photometric properties of the lunar regolith, leveraging the normal albedo map derived from LOLA zero-phase active reflectance data by Lemelin et al. (2016).In this work, we develop an updated calibration procedure and study the data collected through 2021, roughly 8 times more data, including observations of the southern hemisphere, allowing a more detailed examination both globally and of individual targets.This near-infrared data set is unique because of its narrow bandwidth and precise geolocation, and is complementary to observations made by other instruments on board LRO at other wavelengths, including the Diviner Lunar Radiometer Experiment (DLRE; thermal infrared), the Lunar Reconnaissance Orbiter Camera (LROC; visible-ultraviolet), and the Lyman-Alpha Mapping Project (LAMP; far-ultraviolet).

Data Description
LOLA's passive radiometry mode was not in the original instrument mission goals, but was developed in-flight as an alternative mode for when laser altimetric ranging is not feasible (Barker et al. 2016; see Sun et al. 2006 for a similar application of the Mars Orbiter Laser Altimeter).The output of each of LOLA's five detectors is routed through a discriminator designed to detect the return laser pulse in altimetry mode, with all other photon counts recorded as background noise in LOLA's housekeeping data.In its passive radiometry mode, LOLA's laser does not operate, but the sensitivity of the detectors is increased by lowering the discriminator threshold from its normal altimetry setting.In this mode, the solar photons within a bandpass of 1064.45 ± 0.35nm that are reflected off the lunar surface and that are "noise" for altimetry become the signal.This allows LOLA to act as a four-channel radiometer collecting data at 28 Hz (yielding a ∼57 m alongtrack measurement separation for each channel) with a combined footprint diameter of ≈120 m (≈40 m for each channel) from 100 km altitude.The channel corresponding to the center of the footprint, channel 1 of 5, is not used for passive radiometry because it is connected to the laser ranging telescope mounted on LRO's high-gain antenna and thus can also receive photons from the Earth and Sun (Mao et al. 2017).Because LOLA was not designed as a passive radiometer, a calibration procedure (Section 3) is required to obtain absolute radiance from the raw photon counts.
In this work, we analyze results for data collected from 2013 December through 2021 December, over 3.8 billion observations with nearly pole-to-pole coverage and ≈0.01°mean spacing between ground tracks for the northern hemisphere (≈300 m at the equator).Because the first 30 months of data were collected primarily in the northern hemisphere, the mean spacing is ≈30% larger for the southern hemisphere.

Calibration
Before they can be used for scientific analysis, the raw measurements must be calibrated by removing instrumental effects and transforming the receiver count rate to absolute radiance units.Throughout this text, counts refers to the counts measured in a full-rate minor frame, or 1/28th s.
Even without incident light, thermal electrons in the detectors ("dark current") cause nonzero receiver counts.The dark current is a function of detector temperature, plus a small increasing trend over time (≈3% per yr, averaged over the full observation period).Detector temperatures are strongly influenced by reflected radiation from the lunar surface, which depends on LRO's orbital geometry.In order to cover the full temperature range while minimizing the effect of the time trend, we divide our data set's time range into 17 periods (Table A1 in the Appendix), each of which spans the full range of β angle (the angle between the spacecraft angular momentum vector and the Moon-Sun vector).
Over each period, we parameterize the dark current by fitting the nighttime noise counts for each individual channel with a cubic function of temperature T (y = a + bT + cT 3 ; the T 2 term is omitted in order to keep the fit monotonic when extrapolating outside the nighttime temperature range).Figure 1 shows an example, and the full set of parameters is given in Table A2 in the Appendix.The root-mean-square error (RMSE) for the fit is ≈60-70 counts, comparable to the standard deviation of the raw counts at each temperature for the corresponding channel and period.To correct the photometry for the dark current, we subtract the best-fit model for each channel from the raw fullrate daytime noise counts, which can reach >10,000.
The next step is to transform the corrected counts to an absolute radiance scale.This step requires the viewing geometry: The incidence angle i is the angle between the surface normal and surface-to-Sun vectors, the emission angle e is the angle between the surface normal and surface-to-observer vectors, and the phase angle α is the angle between the surfaceto-Sun and surface-to-observer vectors.We calculate these angles using the EMBREE ray-tracing package (Wald et al. 2014),5 the SPICE toolkit (Acton 1996; Acton et al. 2017), 6  and the LOLA 128 pixel-per-degree (ppd) shape model.
For the reference absolute radiance scale, we adopt the nearglobal mosaics from the Kaguya (SELENE) Multiband Imager (MI; Ohtake et al. 2013). 7These mosaics used the photometric correction derived by Yokota et al. (2011) to normalize the MI images to a standard geometry.Thus, for each daytime LOLA full-rate data point, we invert that photometric correction to adjust the MI mosaic's standardized radiance factor (RADF, defined in Section 4) to the LOLA viewing geometry.This is done for the 1049 nm band, as this is the closest MI band to LOLA's 1064 nm wavelength.Although the ratio of adjusted Kaguya RADF/LOLA counts is nearly constant (Figure 2), the responsivity of the LOLA detectors to incident light varies with detector temperature.Using the same time subsets as for the dark-current models, we fit the ratio of Kaguya RADF/LOLA counts for each channel as a linear function of temperature.To obtain a photometrically uniform subset for each fit, we restrict observations to the far side (longitude 90°-270°), midlatitudes (15°-55°), moderate phase angles (15°-55°), and the "medium" albedo group of Yokota et al. (2011).Parameters for the fits are shown in Table A3 in the Appendix; the RMSE of the fits to this restricted data set is ≈1%-2%.We then derive a normalization factor for each channel to ensure consistency between the LOLA passive and active radiometry, which will be used later to study the phase function.This is achieved by computing the median ratio of RADF for the >8000 data points with α < 0.25°to the actively derived LOLA normal albedo obtained by interpolating on the 10 ppd map of Lemelin et al.  (2016). 8This normalization is <6% for all channels (Table A4).We multiply the dark-current-corrected counts for each channel by the calculated RADF/counts ratio, divide by the normalization factor, and take the mean over all four channels to obtain the final calibrated RADF value at each fullrate data point.We use this four-channel full-rate mean RADF to analyze the phase function in subsequent sections.Conservatively estimating the accuracies of the dark current and RADF/counts fits as 2% and 4%, respectively, we estimate the overall accuracy of the full calibration procedure to be ≈5%.The typical uncertainty of the full-rate four-channel mean RADF ranges from ≈1.5% at α = 0°, to ≈2.2% at α = 45°, and to ≈20% at α = 90°.
Figure 3 shows the results of the full calibration for one orbit.Following the dark-current correction, the median nighttime noise counts for each channel have been reduced to absolute values of <20 counts.In Figure 3, high-frequency variations in calibrated RADF are superimposed on a lowfrequency trend of increasing RADF with time of orbit due to the decreasing average incidence angle as the spacecraft approaches the equator.The high-frequency variations compare well with the normal albedo at each data point interpolated from the actively derived 1064 nm LOLA normal albedo map of Lemelin et al. (2016).As a check of the calibration's fidelity, we constructed a model-independent albedo map from the calibrated RADF data set.To compensate for differences in viewing geometry, we normalized each full-rate RADF data point to the median value for its 1°× 1°× 1°three-dimensional pixel, or voxel, of incidence, emission, and phase angle, while restricting to RADF >0.005 and latitudes below 70°to avoid shadows.We regressed this ratio against the LOLA 1064 nm normal albedo map at 4 ppd, obtaining a slope of 0.260 and y-intercept of 0.042 with a mean and standard deviation of residuals of 0.000 and 0.004, respectively.These values loosely correspond to a global mean albedo and the offset between the zero-points of the active and passive photometric data (Barker et al. 2016).Even with new data in the southern hemisphere, our fit is consistent with the values found by Barker et al. (2016, slope 0.27 and intercept 0.04).

Theory
The radiance factor (RADF, also called I/F) is defined as the reflectance relative to a perfectly diffuse surface illuminated vertically: where I 0 (i, e, α) is the observed radiance at 1064 nm at a Sun-Moon distance of 1 au, J 0 = 647 W m −2 μm −1 is the solar irradiance at 1064 nm and 1 au (Thuillier et al. 2003), and D is the Sun-Moon distance in au.RADF can be expressed as where the normal albedo A N is measured at α = 0°and the photometric function F(i, e, α) describes all effects of illumination and viewing geometry (Shkuratov et al. 2011).The effect of A N is clear in our data (Figure 4), which display significant differences in RADF between the brighter highlands and darker maria.Dividing the RADF by the normal albedo removes much (but not all) of the difference between highlands and maria.The result (Figure 5) is called the photometric function, which can be expressed as , , , 3 where the disk function D(i, e, α) describes the brightness distribution over the lunar disk at a given phase angle, and the phase function f (α) describes the overall drop in radiance as α increases.
In this study, we use the Akimov disk function (Akimov et al. 1999;Shkuratov et al. 2011): where the photometric latitude θ and longitude Λ are related to (i, e, α) by ( ) q a = -L i cos cos cos and q = L e cos cos cos , and the parameter ν = 0.34 for maria and 0.52 for highlands.We use these values for the maria and highlands separately unless otherwise noted.Dividing the photometric function by the disk function removes the effect of i and e, leaving Ψ(α) ≡ RADF/A N /D as the data (Figure 6) to be described by the phase function, f (α).
Generally, there are two main approaches to modeling the phase function.One approach (e.g., Hapke 2012a) is based on a detailed solution of the radiative transfer equation applied to light scattering from a particulate surface.One advantage of this method is that the model parameters can correspond, at least in principle, to specific physical processes or characteristics.In practice, however, this correspondence is not always clear, and finding a unique solution can be difficult when fitting observations (Shepard & Helfenstein 2007;Sato et al. 2014).This could be due, in part, to strong correlations between the parameters and the high number of parameters (up to nine, though this can be reduced with simplifying assumptions; Hapke 2012a).The second approach is more empirical and uses a simplified theory with fewer, less physically specific parameters, but this type of phase function provides relatively unambiguous fits to data in practice (Shkuratov et al. 2011).
Here, we will take the latter approach with the new data set and leave the former for future work.
The first phase function that we use was proposed by Korokhin et al. (2016): where η gives the slope and ρ the "bend" of the phase curve.Velichko et al. (2022) used this phase function to analyze the effects of the landing jets on the surface at the Luna 23 and 24   4)).The increasing spread in values at higher phase angles is caused by small errors in the calculated photometric angles, with the concentration of low values resulting from topographic shadows.
landing sites, finding lower values of ρ in disturbed areas, which they attributed to destruction of the porous "fairy-castle" microstructure and surface smoothing due to deposition of dust.
We subsequently refer to Equation (5) as the two-parameter phase function.
The second phase function that we will use is described in Shkuratov et al. (2011): 1 2 where m and μ 1 are inversely correlated with the relative amplitude and angular width, respectively, of the opposition surge (increasing radiance at very low α), and μ 2 (which is hypothesized to be correlated with surface photometric roughness) gives the overall slope of the phase curve.We note that photometric roughness is thought to be dominated by roughness at the smallest scale at which shadowing exists, typically on the order of 0.1 mm for the Moon (Shepard & Campbell 1998;Helfenstein & Shepard 1999).To avoid ambiguity in the fitting procedure, we enforce μ 1 > μ 2 .Using this phase function with Earth-based optical telescope observations, Shkuratov et al. (2011) found good correspondence between the parameters and lunar geological units.We subsequently refer to Equation (6) as the three-parameter phase function.
To assist with intuition, we illustrate the effects of varying the parameters of both phase functions in Figure 7.For the twoparameter phase function (Equation ( 5)), we see that increasing the slope η (while holding ρ constant) lowers the phase function at all α by an amount that increases with α (Figure 7(a)), while increasing the bend (decreasing ρ while holding η constant) steepens the phase function at very low phase angles (α 5°) while flattening it at all other α (Figure 7b).For the three-parameter phase function (Equation ( 6)), there is a clear distinction between the ranges of phase angles affected by μ 1 and μ 2 , while m affects the overall shape by weighting the low-and high-α components (Figures 7(c)-(e)).Although Equations ( 5) and (6) are not equivalent, we will see in Section 4.2.2 that ρ roughly corresponds to the opposition effect inverse amplitude m, but  Notes.For the two-parameter phase function (Equation ( 5)), η gives the overall slope and ρ gives the bend.For the three-parameter phase function (Equation ( 6)), μ 1 gives the width of the opposition effect, μ 2 gives the overall slope, and m gives the relative amplitude of the opposition effect.Error bars are ±2× standard error of the fit, using the standard error of the mean as the measure of uncertainty in the data.
is unrelated to the opposition effect width μ 1 .Only the threeparameter phase function is able to model the opposition effect width at low α without affecting values at higher α (Figure 7(c)).

Fitting the Phase Function
For all phase-function fits, we restrict the data to i < 70°, e <60°, and Ψ(α) > 0.05 to minimize the effect of shadows.We build a histogram of the full-rate data points in 1°× 0.01 bins of α and Ψ(α), and use the median value of Ψ(α) for each α as the data to be fitted.We carry out the fits with the lmfit Python package for nonlinear least-squares minimization,9 using robust regression with Huber weighting (Huber 1981).

Global Data
First, we compute fits on a global scale, while treating the highlands and maria separately (Figure 8, Table 1) using the maria boundaries of Nelson et al. (2014).To produce the global phase curves, we first build histograms of the full-rate data points for each 1°× 1°map pixel equatorward of 70°and find the median value of Ψ(α) for each α.We then take the median of the pixel medians at each α as the data to be fitted.This procedure effectively normalizes variations in data density over latitude and phase angle, reducing any biases that would result from building a single histogram with all of the full-rate data points.The two-parameter phase function yields a mean residual of 4.7% and standard deviation of 4.6% for latitudes equatorward of ±70°.The three-parameter phase function yields comparable residuals, with mean 4.1% and standard deviation 5.0%.
Overall, the phase function for the highlands drops off more slowly than for the maria (Figure 8).For the two-parameter function, the parameter η has a higher value for the maria, while ρ has a higher value for the highlands.This is consistent with the results of Korokhin et al. (2016), whose analysis of LROC observations of Mare Nubium found η to be inversely correlated with albedo.Comparison with the results of Velichko et al. (2022) suggests that the higher value of ρ in the highlands indicates greater photometric roughness.For the three-parameter function, the values of μ 1 and m indicate a broader and more pronounced opposition effect for the maria, and μ 2 suggests greater photometric roughness for the highlands than for the maria.We further discuss the interpretation of ρ and μ 2 in terms of photometric roughness in Section 7.  6)) with separate parameters for the maria and highlands (Table 1).Next, we construct a photometrically normalized near-global mosaic of the LOLA passive radiometry.We normalize the RADF of each full-rate point to a standard geometry, RADF , , 30 , 0 , 30 30 , , , using the maria or highlands phase function as appropriate.
Figure 9 shows the resulting map at 16 ppd adopting the threeparameter phase function.

Equatorial Strip
As an alternative to fitting the global highlands and maria data, we now perform a spatially resolved analysis of the phase function over a narrow equatorial strip (±3°latitude) because this is where we have the most complete coverage of low phase angles (Figure 10), which are particularly important for obtaining a reliable three-parameter fit.For example, in the map pixel shown in Figure 11, the differences in parameter values when the α < 20°data are excluded in the threeparameter fit are significant, especially for μ 1 , and the restricted  fit is noticeably worse for α < 20°.This is expected since the purpose of μ 1 is to describe the opposition surge at low α.
Fitting the two-and three-parameter phase functions for each pixel in the equatorial strip (Figure 12) confirms previous studies that found systematic trends of the phase function with albedo (Yokota et al. 2011;Sato et al. 2014;Barker et al. 2016).Although the measured phase function Ψ(α) is normalized by normal albedo, it remains sensitive to other regolith characteristics correlated with albedo, as we discuss later.
Next, we fit each phase-function parameter versus albedo using a continuous piecewise linear function with segments corresponding to typical maria and highlands albedo values (Figure 13).We then evaluate the photometric function at all full-rate data points globally, using these piecewise fits to generate the phase-function parameters from the normal albedo at each point.Calculating residuals globally, we find that the albedo-based phase-function fits (Table 2) are more closely centered on zero and have lower standard deviations (RMSE of 0.3% ± 3.3% for two parameters and 0.7% ± 3.4% for three parameters) compared to the globally derived phase-function fits in Table 1.This indicates that the phase function averaged over ∼30 km scales can be simplified and made more accurate if parameterized by the LOLA normal albedo of Lemelin et al. (2016).
Comparing the phase-function parameter fits for each map pixel also provides insight regarding the relationship between the two-and three-parameter phase functions.The strongest correlations are between m and both η (correlation coefficient, R = −0.90)and ρ (R = 0.91), all three of which are very strongly correlated with A N (R = 0.92, −0.96, and 0.92 for m, η, and ρ, respectively; see Figure 13).This suggests that these parameters are associated with surface properties that correlate with albedo, such as composition and exposure age (discussed further in Sections 5.1 and 7).For the maria pixels (A N < 0.22), the roughness parameter μ 2 , which dominates the overall slope of the three-parameter phase function as α increases, is moderately correlated with the two-parameter slope η (R = 0.64) and normal albedo (R = −0.47).However, for the highlands pixels (A N > 0.22), μ 2 has only a weak correlation with η (R = 0.30) and is uncorrelated with A N (R = 0.03).The opposition surge width μ 1 is not correlated with η or ρ (R = −0.03,0.02) for either highlands or maria.There is also no correlation (R = 0.03) between μ 1 and A N , agreeing with the theory that normal albedo is affected by the amplitude but not the width of the opposition surge (Hapke 2012a).This suggests that μ 1 is largely controlled by, and μ 2 is at least influenced by, surface properties unrelated to albedo.We discuss the apparently significant effect of A N on the phase function in Section 7, along with physical interpretations of the phasefunction parameters.Global maps of the A N -dependent parameters can be derived by applying the fits in Table 2 to the normal albedo map.

Model-independent Analysis
Several geologically important surface characteristics of the optically active portion of the regolith can influence the phase function at 1064 nm (McGuire & Hapke 1995;Barker et al. 2016).Here, we consider FeO, optical maturity (OMAT; Lucey et al. 2000b), submicroscopic iron (SMFe; Trang & Lucey 2019), surface slope, and TiO 2 (Lucey et al. 2000a).Consistent with Barker et al. (2016), we find that LOLAderived 50 m baseline roughness (Neumann et al. 2015) and rock abundance (Bandfield et al. 2017) have little effect on Ψ (α) when controlled for OMAT.We consider Diviner-derived infrared thermophysical properties later, in the discussion of the Chaplygin cold spot (Section 6.2).Here, we restrict our analysis to within ±50°latitude because the Kaguya MI maps of FeO and OMAT (Lemelin et al. 2019) cover only that range.
In Figure 14, for both highlands and maria, increasing FeO is associated with decreasing values of the phase function (the difference at α = 30°between the plotted bins containing 99% of the data, defined as Δ 30 for brevity, is −0.09 for both highlands and maria).Although the FeO distribution is concentrated at lower values for highlands and higher values for maria, the phase functions are very similar (the difference at α = 30°between the corresponding highlands and maria bins is <0.005 for 5 <FeO <17 wt%).FeO has a strong negative correlation with the 1064 nm normal albedo (R = −0.95),consistent with the well-known absorption minima of minerals containing ferrous iron (Fe 2+ ) at around 1 μm (e.g., Horgan et al. 2014).
In Figure 15, the phase-function values increase with OMAT, with a markedly larger effect for the highlands than the maria (Δ 30 = 0.11 and 0.03, respectively).Since the OMAT parameter decreases with maturity, this represents a decrease in the phase function as space weathering progresses.We also note that OMAT is designed to be independent of FeO content, although in practice there is a small correlation, so that OMAT is best compared between areas of similar FeO (Lucey et al. 2000b).It is likely that the apparently larger effect of OMAT on the highlands phase function contains some contribution from FeO, because the least mature highlands tend also to have low FeO.Restricting the highlands data to FeO >5 wt% does reduce the effect of OMAT, but the difference in phase-function values remains significant (Δ 30 = 0.08).
Exposure to the space environment has a number of significant effects on an airless surface, including changes in particle size distribution, the production of glassy agglutinates, and the accumulation of SMFe on the rims of mineral grains and within agglutinates (Hapke 2001;Trang & Lucey 2019).Given that the phase function is affected by both FeO (the source of SMFe particles) and OMAT, we would expect it to be affected by SMFe as well.In Figure 16, the phase-function values decrease with SMFe, with a larger effect for the highlands than the maria (Δ 30 = −0.09and −0.06, respectively).However, because SMFe depends on FeO and is correlated with OMAT (Trang & Lucey 2019), this is not entirely independent information.The effect of SMFe largely disappears if we limit FeO and OMAT to the peaks of their distributions (FeO: 5-7 wt% highlands, 16-18 wt% maria; OMAT: 0.14-0.18)Conversely, if we limit SMFe to the peak of its distribution (1.3-1.5 wt%), we find that both FeO and Notes.For the two-parameter phase function (Equation ( 5)), η gives the overall slope and ρ gives the bend.For the three-parameter phase function (Equation ( 6)), μ 1 gives the width of the opposition effect, μ 2 gives the overall slope, and m gives the relative amplitude of the opposition effect.
OMAT still significantly affect the phase function when the other is also restricted.This suggests that SMFe captures some, but not all, of the effect of optical maturity on the 1064 nm phase function.
Phase-function values increase similarly for highlands and maria with topographic surface slope (Δ 30 = 0.06 for both; see Figure 17), although slopes are generally much lower for the maria.This behavior may result from mass wasting on crater walls exposing optically immature material (Barker et al. 2016), regardless of the regolith FeO content.Consistent with this hypothesis, restricting OMAT to the peak of its distribution (≈0.14-0.18)reduces Δ 30 to 0.02 for the highlands.
There is a moderate effect of TiO 2 for maria (Δ 30 = −0.07;see Figure 18), but the highlands phase function is very similar to the maria at >1-2 wt%.The small percentage of highlands data points with higher TiO 2 primarily lie near maria boundaries, where terrain may be misclassified and where mixing with mare material due to impact ejecta is probable (Li & Mustard 2005).In the maria, the effect of TiO 2 saturates at 8 wt%.

Partial Least-squares Analysis
We seek to quantitatively measure which geological properties are the dominant controls on the phase function.Determining the relative contribution of each property to the phase function requires a regression method that can account for correlations among the independent variables.(For simplicity, we assume that the relationships are at least approximately linear.)Partial leastsquares regression (PLS; Abdi 2003) is similar to principal component analysis (PCA) in that both methods decompose the problem into orthogonal components; however, while PCA chooses components to best explain the independent variables (e.g., FeO, OMAT), PLS chooses components to best explain the covariance between the independent and dependent variables (e.g., ρ, η).
For our analysis, we bin the phase-function data into threedimensional voxels of FeO, OMAT, and either surface slope (highlands) or TiO 2 (maria), dividing each independent variable into 10 bins based on the 1st to 99th percentiles of its global values.Choosing this number of bins and variables is intended to include the geological properties with the largest effect while having enough data points per voxel to properly resolve the phase function.We omit slope for maria and TiO 2 for highlands because of their small range of values; see Figures 17 and 18.In order to avoid outliers and misclassified points, we consider only the FeO bins that contain >1% of the data (3.40-13.48wt% and 8.44-20.20 wt% for highlands and maria, respectively; see Figure 14).We then fit the two-parameter phase function for each voxel, and take η and ρ to be the dependent variables for PLS.We found that the high sensitivity of the three-parameter phase function to the absence of low-α data (Figure 11) prevented a reliable PLS analysis.Because not all voxels contain a broad range of phase angles, which may lead to inaccurate fits, we keep only voxels for which the relative uncertainty of both parameters is <2% and the RMSE of the fit is <0.03.Before performing the regression, we standardize each variable to have zero mean and a standard deviation of 1, which makes it easier to compare their relative contributions.The results are shown in Figure 19 and Table 3.For both highlands and maria, the variance of η is well explained (≈98% for highlands and ≈93% for maria), with the largest effect from OMAT in the highlands and from FeO in the maria (as given by their larger absolute standardized coefficients) and the smallest from slope/TiO 2 .The variance of ρ is not as well explained (only ≈87% for highlands and ≈63% for maria), and is dominated by FeO for both highlands and maria.In the maria, the effect of TiO 2 on the phase function (and thus on the twoparameter fit) decreases at higher TiO 2 (Figure 18), and this departure from linearity may help account for the lower variance of ρ explained by PLS.These results confirm that the individual effects of FeO, OMAT, surface slope, and TiO 2 (Figures 14-15 and Figures 17-18) remain significant when multiple geological properties are considered simultaneously.

Reiner Gamma Swirl
The Reiner Gamma swirl is a photometric anomaly located in western Oceanus Procellarum, characterized by increased albedo (Figure 20) and unusual spectral properties across a broad range of wavelengths, and coincident with a crustal magnetic anomaly (Hood & Schubert 1980).On-swirl areas have been found to have spectral characteristics similar to fresh impact craters and ejecta (Denevi et al. 2016) but inconsistent with basaltic immature regolith (Pieters et al. 2021).There are several leading hypotheses regarding swirl formation, including that the space weathering of the swirl material may be unusual due to the magnetic anomaly deflecting the solar wind, but not micrometeoroids (e.g., McFadden et al. 2019).Under this "solar wind stand-off" hypothesis, swirl material likely has lower concentrations of SMFe, although studies disagree about the affected size ranges of SMFe particles (Trang & Lucey 2019; We compare the phase function between the brighter onswirl areas and the darker surrounding off-swirl maria, finding distinctly higher values for the on-swirl data at all phase angles (Figure 21).The two-parameter phase function shows a lower slope (η) for the on-swirl region, consistent with its higher albedo, while the difference in the bend (ρ) of the phase function is negligible.The three-parameter phase function shows a lower-amplitude opposition surge (i.e., higher m) with similar width (difference in μ 1 < standard error of fit) and lower roughness (smaller μ 2 ) for the on-swirl regions (Figure 21, Table 4).
Possible interpretations of the phase-function parameter differences depend on whether the values of FeO and OMAT are accurate for swirl material, which may experience an atypical mixture of space-weathering processes.The on-swirl areas have lower FeO (14.7 ± 1.0 versus 17.7 ± 0.8), higher OMAT (0.25 ± 0.03 versus 0.18 ± 0.02), and lower TiO 2 (3.8 ± 0.9 versus 6.0 ± 1.1) than the surrounding maria, consistent with the higher phase curve of the on-swirl areas (Section 5.1).However, the FeO/TiO 2 /OMAT algorithms derived by Lucey et al. (2000b) and used by Lemelin et al. (2019) relied on Apollo samples, which were exposed to both solar wind and micrometeoroid impacts, so the solar wind stand-off hypothesis suggests that their calibration may not be applicable to the spectrally unusual swirl material.
We first consider the possibility that FeO, OMAT, and TiO 2 are unreliable on swirls, and interpret the phase-function parameters assuming that the on-swirl material is compositionally similar to the surrounding maria.In this case, the   If instead we take the differences in FeO, OMAT, and TiO 2 between on-and off-swirl areas at face value, our data show that the phase curves for both on-and off-swirl areas align with those of compositionally equivalent ("on-swirl-like" and "offswirl-like") global maria (Figure 22).This differs from Barker et al. (2016), who found that the on-swirl phase-function values for the first year of data were higher than those of its compositionally equivalent maria, as defined using Clementine FeO, OMAT, and TiO 2 .However, Clementine FeO is ≈2 wt% higher than Kaguya FeO on-swirl but comparable off-swirl, and the relatively small sample of compositionally equivalent northern-hemisphere maria in the earlier study may not be representative globally.Indeed, the on-swirl phase-function values for the first year of data do align with the global onswirl-like phase curve for the full 8 yr.
Figure 22 suggests that swirl material has a similar spaceweathering history to that of comparable maria outside of swirls, which disagrees with the solar wind stand-off hypothesis.As we discuss in Section 7, multiple scattering reduces the effect of photometric roughness, so a different interpretation of our results that does not depend on surface texture is that the lower on-swirl value of μ 2 is influenced by stronger multiple scattering due to its higher albedo.Similarly, the higher onswirl value of m is consistent with the higher albedo of the swirl material because the amplitude of the opposition effect is inversely correlated with single-scattering albedo (Sato et al. 2014).This interpretation agrees with Kinczyk et al. (2022), who used LROC Narrow Angle Camera (NAC) imagery and the Hapke (2012a) model to find that differences in photometric behavior between on-and off-swirl regolith could be explained by the higher single-scattering albedo of the swirl material, without any variations in porosity or roughness.Domingue et al. (2023), also using LROC-NAC imagery and the Hapke (2012a) model, likewise found no differences in microscale roughness between on-and off-swirl areas in the Mare Ingenii swirl region.While it remains possible that there are also textural differences between the on-and off-swirl surfaces, the data presented in Figure 22 suggest that the phase curves depend primarily on composition and optical maturity.

Chaplygin Cold Spot
Diviner observations have been used to identify "cold spots" surrounding recently formed craters distributed over the lunar surface (Bandfield et al. 2014;Williams et al. 2018) using the latitude-adjusted regolith nighttime temperature (ΔT; Powell et al. 2023) and the H parameter (a measure of near-surface thermal inertia and porosity; Hayne et al. 2017).One of the largest of these cold spots is located near Chaplygin crater, covering a mostly Nectarian-aged area of the equatorial highlands.We hypothesize that the porosity variations inferred in the top ≈10 cm of regolith at the cold spots may be observable in the 1064 nm phase function.To test this hypothesis, we delineate this cold spot using ΔT, and compare its phase function with that in a surrounding annulus that does not contain any cold spot material (Figure 23).The two areas are similar in geological surface characteristics other than ΔT and H, such as FeO, TiO2, OMAT, normal albedo, and surface slope.We find little to no difference in the phase function between the cold spot and its surroundings (Figure 24), even though both ΔT and H appear to influence the global phase function when considered alone.However, both ΔT and H are globally correlated with OMAT (R ≈ 0.5 and −0.4, Figure 20.The Reiner Gamma swirl seen in normalized RADF from LOLA passive radiometry (32 ppd).We define on-swirl as data points lying within both the yellow bounding box and the blue swirl outlines (shapefiles from Denevi et al. 2016), and off-swirl as data points lying inside the area plotted here and outside both the on-swirl bounding box and the swirl outlines.
respectively), and when we restrict OMAT to the range seen at Chaplygin, the ΔT and H influence on the global phase function disappears.We observe similar results for other cold spots, and thus conclude that the cold spots detected by Diviner are not visible in the 1064 nm phase function, possibly due to the difference between LOLA's sensing depth (likely 1 mm), and the ≈10 cm depth over which porosity is inferred from Diviner temperature data.This result is consistent with the lack of a clear signature of cold spots in maps of albedo from Clementine ultraviolet/visible/near-infrared and LROC-NAC images (Bandfield et al. 2014).Speyerer et al. (2016) found that the top 2 cm of regolith could be fully overturned by secondary cratering in 81,000 yr, which suggests that the surface microstructure could return to equilibrium over a much shorter  Notes.For the two-parameter phase function (Equation ( 5)), η gives the overall slope and ρ gives the bend.For the three-parameter phase function (Equation ( 6)), μ 1 gives the width of the opposition effect, μ 2 gives the overall slope, and m gives the relative amplitude of the opposition effect.Error bars are ±2× standard error of the fit, using the standard error of the mean as the measure of uncertainty in the data.
timescale than the 100 kyr to 1 Myr persistence time of a typical cold spot (Williams et al. 2018).As cold spot formation may not involve significant deposition or erosion of surface material (Bandfield et al. 2014), restoration of the surface microstructure would remove any difference in properties known to affect the 1064 nm phase function, rendering cold spots undetectable by LOLA.

Silicic Anomalies
The Compton-Belkovich Volcanic Complex (CBVC) is an area of silicic volcanism in the northern highlands of the lunar far side (60.5°N, 99.5°E) that is associated with a high-thorium anomaly (Jolliff et al. 2011).CBVC consists of an irregularly shaped central caldera surrounded by volcanic domes to the north, east, and west.Using data from the Mini-RF S-band radar on board LRO, Petro et al. (2013) found that much of CBVC is covered by a mantle of fine-grained pyroclastic material.Chauhan et al. (2015) classified CBVC as an ash-flow caldera characterized by late-stage pyroclastic eruptive phases.Lucey et al. (2000b) found generally low OMAT values for basaltic pyroclastic deposits on the central nearside, attributing this not to a high exposure age but instead to glassy material that violates the assumptions underlying their OMAT algorithm.In contrast, CBVC stands out from its surrounding highlands as a high-OMAT (low exposure age) anomaly.Since the eruption style at CBVC may have been different from other silicic anomalies (Head & Wilson 2017), it is possible that its phase function is different, too.
In Figure 25, we compare the CBVC phase function with those for several silicic anomalies at lower latitudes that are embayed by maria and characterized by an extrusive eruption style, including Darney-τ (−11.84°,333.32°),Gruithuisen-γ (36.56°, 319.28°),Gruithuisen-δ (36.07°, 320.41°), and the Helmet (−16.61°,328.88°).Clegg-Watkins et al. (2017) used LROC-NAC imagery and the Hapke (2012a) model to infer higher reflectance and higher single-scattering albedo at CBVC than at the Gruithuisen domes, suggesting a lower mafic mineral content and the presence of silicic glassy beads in the CBVC pyroclastic material.We observe much higher 1064 nm phase-function values at CBVC than at the other silicic anomalies (Figure 25), with lower η, consistent with CBVC's higher normal albedo.At high α, the phase functions are consistent with the FeO/OMAT-equivalent global highlands data.Note that the high latitude of CBVC places it outside the coverage of the Kaguya FeO and OMAT maps, so we use the corresponding Clementine maps at these sites (Lucey et al. 2000a).However, a strip of missing data in the Clementine maps means that we have FeO and OMAT values for only the western half of CBVC, although the LOLA ground tracks are spread across the full area.
At α < 60°, the observed phase-function values for CBVC are generally below those of compositionally equivalent global highlands (Figure 25).Due to the high latitude of CBVC, observations at α < 60°required an LRO targeting slew, and these ground tracks mostly cover the central caldera region for which the Clementine maps have no data.One exception is a track at α ≈15°, which crosses the northern dome and has a phase-function value consistent with the compositionally equivalent global highlands.These differences in the 1064 nm phase function suggest a different FeO and/or OMAT between the central caldera and the west dome.From the data plotted in Figures 14 and 15, we estimate that the ≈0.04-0.05decrease in the CBVC phase function (solid points) relative to its compositionally equivalent highlands (solid line) in Figure 25 for 20°< α <55°is approximately consistent with   (Powell et al. 2023).We define our cold spot sample as the ∼167,000 LOLA data points inside the inner yellow circle for which ΔT is below the global 1st percentile (−3.65°C, assigned a color of white).The control sample is the ∼222,000 LOLA data points between the yellow circles.an increase of 3-4 wt% FeO or decrease of 0.06-0.08 in OMAT if either property is considered alone.
An alternative, and perhaps more likely, interpretation is that the disk function (Equation ( 4)) is less accurate for the slewed data due to their higher emission angle compared to the nadir data (median e = 44.0°versus2.7°).Indeed, we find much better agreement between the phase function of the CBVC slewed data (α < 60°) and that of compositionally equivalent highlands when restricting the latter to similarly high emission and off-nadir angles.This supports the interpretation that FeO and OMAT are the primary controls on the highlands phase function for silicic anomalies.

Discussion
Both the 1 μm absorption band associated with minerals containing ferrous iron (e.g., olivine, clinopyroxene, orthopyroxene) and the darkening caused by space weathering affect the 1064 nm albedo, but we find that normalizing RADF by the normal albedo does not eliminate the effects of FeO and OMAT on the phase function.This behavior suggests that surface properties correlated with normal albedo (such as multiple scattering and the backscattering behavior of individual particles) affect the phase function.There is some ambiguity associated with the semiempirical phase functions we use here because their parameters can be interpreted in terms of both surface microtexture and scattering properties of individual particles.However, as we discuss next, the correlations with normal albedo found in Section 4.2.2 suggest that η, ρ, and m are dominated by individual particle properties and that μ 1 and μ 2 are dominated by texture.
High photometric surface roughness, i.e., the mean height of surface irregularities relative to their mean spacing, reduces reflectance at mid-to-high α by increasing shadowing (Korokhin et al. 2007).Theoretical considerations indicate that increased porosity may narrow the opposition peak, as does an increase in the width of the grain size distribution (Helfenstein & Shepard 2011;Hapke 2012a).Roughness and porosity variations may be localized and short lived, however, due to the ubiquity and short timescales of micrometeoroid bombardment and small secondaries (Sato et al. 2014;Speyerer et al. 2016).Indeed, Sato et al. (2014) produced good fits to observations using the Hapke (2012a) model while holding the photometric roughness constant over large areas of the lunar surface.Therefore, μ 1 and μ 2 would be expected to be approximately constant for most of the surface, consistent with what we find.Because of the porous microstructure of the regolith surface, multiple scattering can illuminate internal shadows and increase reflectance across a broad range of α, countering the effect of surface roughness.The strength of multiple scattering increases with albedo, which is consistent with the correlations we find between normal albedo and η, ρ, and m.However, the weak correlation between μ 2 and normal albedo suggests that photometric roughness is the dominant controlling factor for this parameter.
Increasing FeO, increasing TiO 2 , and decreasing OMAT all decrease the phase curves across a broad range of α.This may be due to changes in the relative strength of forward-and backscattering by individual grains (i.e., the single-particle phase function), which depends on grain properties such as composition, size, surface roughness, internal scattering, and internal absorption (McGuire & Hapke 1995;Sato et al. 2014).Interpolating LOLA data points onto the mineral maps of Lemelin et al. (2019), we see that the downward trend in the phase function with increasing iron abundance corresponds to increasing concentrations of the mafic minerals olivine, orthopyroxene, and clinopyroxene; conversely, the phase function shows an upward trend as the concentration of plagioclase feldspar increases.The smaller downward trend in the phase function with increasing TiO 2 in the maria is associated with higher concentrations of the opaque mineral ilmenite, the dominant carrier of titanium in lunar materials (Lucey et al. 2000a).With the scattering behavior of single particles described by a double-lobed Henyey-Greenstein phase function (McGuire & Hapke 1995), and using the "hockey-stick relation" (Hapke 2012b) to determine the balance between backscattering and forward scattering, the differences in the phase curves are qualitatively consistent with a decrease in backscattering with increasing FeO, exposure age, and TiO 2 .This would agree with the results of Sato et al. (2014), who found that backscattering decreases with increasing olivine, pyroxene, ilmenite, SMFe, and/or low-albedo agglutinates, although testing this would require analysis of our data using the Hapke (2012a) model.Both the 1 μm absorption minima of mafic minerals and the opacity of ilmenite result from the low albedo of individual particles, caused by internal absorption for the mafic minerals and low Fresnel reflection at grain surfaces for ilmenite (Hapke 2001).This suggests that as the single-particle albedo decreases (e.g., FeO increases), the backscattering strength decreases, which may contribute to the correlations we find between normal albedo and η, ρ, and m.The Hapke model's single-scattering phase function may contribute at some level to the overall slope, μ 2 .However, as explained above, the observed weak dependence of μ 2 on normal albedo suggests photometric roughness has a larger effect on μ 2 than does the single-scattering phase function.The weakening of multiple scattering with decreasing singleparticle albedo could also be a contributing factor to the observed parameter correlations with normal albedo, and further analysis will be needed to differentiate the effects of single and multiple scattering.
Finally, the opposition surge width increases with decreasing OMAT (higher exposure age) in the highlands and immature maria.This would be consistent with lower porosity and/or a narrower distribution of grain sizes as the grains of fresh ejecta are broken down by micrometeoroid impacts before equilibrium with agglutination is reached (McKay et al. 1991;Sato et al. 2014).However, the opposition surge in the mature maria narrows with decreasing OMAT.Sato et al. (2014) suggest that this effect depends on mineral composition, because differing rates of comminution between titanium-containing ilmenite, which is more common in the maria, and other silicates may instead produce a broader grain size distribution as space weathering progresses.
The simplified phase functions we use here fit our data reasonably well, but further investigation of the effects of specific processes will likely require analysis with a more detailed photometric model such as that of Hapke (2012a).

Summary
We have presented a newly calibrated near-infrared passive radiometry data set that expands the data analyzed by Barker et al. (2016) to include approximately 8 times more observations with coverage of both hemispheres.The calibration converted LOLA's passive radiometry to absolute radiance by removing the effects of thermal noise in the detectors and referencing to the Kaguya MI mosaic.We fitted the data to two phase functions based on simplified physics to analyze the global highlands and maria and create a 1064 nm reflectance map at standard viewing geometry.The albedo-dependent phase functions we derived here can be used by other investigators to predict the near-infrared surface reflectance even in the polar regions given the map of LOLA normal albedo (Lemelin et al. 2016).
Examining the effect of individual geological properties, we found that the phase curves of the data are strongly affected by FeO and OMAT, and moderately affected by surface slope and TiO 2 .However, the thermophysical properties observed in the infrared by Diviner appear to have no effect on the 1064 nm phase function, possibly due to the LOLA sensing depth being much shallower than the ≈10 cm depth over which the Diviner properties are inferred.By applying PLS regression to phasefunction fits, we quantitatively confirmed that FeO and OMAT are the dominant influences on the phase function when multiple geological properties are considered together.Examining the Reiner Gamma swirl and several silicic anomalies, we found results generally consistent with varying FeO and OMAT.Further work is needed to examine how other properties that vary with iron abundance and exposure age, such as multiple scattering and the scattering behavior of individual grains, affect the 1064 nm phase function.

Appendix
Tables A1-A4 are for the calibration procedure described in Section 3.

Figure 2 .
Figure 2. Ratio of Kaguya RADF to LOLA counts for Channel 2 from 2021 January 31 (β = −85.8°) to 2021 July 19 (β = 86.9°).The temperature range of the data is smaller than in Figure 1 because the range of latitude (which affects detector temperature) is restricted.The black line shows best-fit linear function of temperature (T) for this time range (RMSE ≈1.7 × 10 −6 , or ≈5.7% of the median ratio for this time range).

Figure 4 .
Figure 4.The histogram of RADF for the entire data set (1°× 0.005°bins) shows the effect of albedo differences between brighter highlands and darker maria, with the concentration of low values with RADF <0.005 resulting from topographic shadows.

Figure 5 .
Figure 5.The photometric function for the entire data set (1°× 0.01°bins) is explicitly corrected for normal albedo and contains all information about illumination and viewing geometry.

Figure 6 .
Figure 6.Phase function for the entire data set (1°× 0.01°bins), taking ν = 0.43 in the disk function (Equation (4)).The increasing spread in values at higher phase angles is caused by small errors in the calculated photometric angles, with the concentration of low values resulting from topographic shadows.

Figure 8 .
Figure 8. Fits to the phase function for global highlands (top) and maria (bottom).Yellow circles show the median values in 1°× 1°map pixels and 1°bins of phase angle.

Figure 10 .
Figure 10.Minimum phase angle in 1°× 1°map pixels.The minimum possible phase angle while nadir-pointing is dependent on latitude, increasing poleward.Stripes of lower minimum phase at mid to high latitudes result from slews in which LRO was pitched ≈45°forward on every other orbit to support LROC imaging.

Figure 11 .
Figure 11.Phase-function fits for the equatorial strip pixel at longitude 27°-28°, latitude 0°-1°(in Mare Tranquilitatis) illustrate the high sensitivity of the threeparameter phase function to a lack of low-α data.The orange and blue lines show the fits to all data for the two-and three-parameter functions, respectively, and the green line shows the fit for the three-parameter function with the data restricted to α > 20°.The two-parameter fit using only the α > 20°data is not shown because the parameter values change so little (η = 1.1985, ρ = 0.4547) that the results plot on top of the all-α fit.

Figure 13 .
Figure 13.Piecewise linear fits between A N and phase-function parameters for the ±3°latitude equatorial strip (Figure 12), with marginalized histograms for each parameter.The break at A N = 0.22 divides typical maria (<0.22) and highlands (>0.22) albedo values, as seen in the marginalized histogram for A N .

Figure 14 .
Figure 14.Phase-function dependence on FeO, with bins ranging from the 1st to 99th percentile of all full-rate data points interpolated on the Kaguya MI FeO map at 128 ppd.Curves show median over 1°bins of phase angle.Percentages in legend indicate proportion of the plotted data in each FeO bin for highlands and maria, respectively.Nonsmooth curves result from scarcity of data in the corresponding bins.

Figure 15 .
Figure 15.Phase-function dependence on OMAT, with bins ranging from the 1st to 99th percentile of all full-rate data point values interpolated on the Kaguya MI OMAT map at 128 ppd.Curves show median over 1°bins of phase angle.Percentages in legend indicate proportion of the plotted data in each OMAT bin for highlands and maria, respectively.

Figure 16 .
Figure 16.Phase-function dependence on SMFe, with bins ranging from the 1st to 99th percentile of all full-rate data point values interpolated on the Trang & Lucey (2019) map at 30 ppd.Curves show median over 1°bins of phase angle.Percentages in legend indicate proportion of the plotted data in each SMFe bin for highlands and maria, respectively.

Figure 17 .
Figure 17.Phase function by surface slope, with bins ranging from the 1st to 99th percentile of all full-rate data point values interpolated onto a map calculated from the LOLA digital elevation model at 128 ppd.Curves show median over 1°bins of phase angle.Percentages in legend indicate proportion of the plotted data in each surface slope bin for highlands and maria, respectively.

Figure 18 .
Figure 18.Phase function by TiO 2 , with bins ranging from the 1st to 99th percentile of all full-rate data point values interpolated onto the Clementine TiO 2 map (Lucey et al. 2000a) at ≈200 m pixel -1 .Curves show median over 1°bins of phase angle.Percentages in legend indicate proportion of the plotted data in each TiO 2 bin for highlands and maria, respectively.

Figure 19 .
Figure 19.PLS fits of the parameter values obtained from fitting the two-parameter phase function for each (FeO, OMAT, slope/TiO 2 ) voxel vs. the corresponding values of the geological properties of interest.PLS coefficients are shown in the units of each variable rather than the standardized values discussed in the main text.
weakening of the opposition surge and decreased roughness in the on-swirl area collectively suggest a smoother and/or more compacted surface because compaction reduces shadowing at scales that affect reflectance at low α (Hess et al. 2023), consistent with the results of some earlier studies.Bhatt et al. (2023) used Earth-based imaging polarimetry to suggest smoothing and compaction of the upper layer of the regolith within the swirl, while Hess et al. (2023) reached a similar conclusion by comparing the swirl to the effects of the landing rockets at the Chang'e-5 site.This interpretation is inconsistent with the electrostatic dust transport hypothesis of Garrick-Bethell et al. (2011), because redeposition of lofted dust is associated with formation of a highly porous "fairy-castle" microstructure (Byron et al. 2019).However, Garrick-Bethell et al. (2011) hypothesized that the transported fine dust is brighter and more feldspathic (lower iron abundance) than the surrounding maria, resulting in a different composition for the swirl material.

Figure 21 .
Figure 21.Fits to the phase function for on-swirl (top) and off-swirl (bottom) areas of Reiner Gamma.Median of data over 1°bins of phase angle is shown by yellow circles.

Figure 22 .
Figure 22.Median (circles) and interquartile range (shaded) of the phase function in the vicinity of the Reiner Gamma swirl (Figure20).The phase function is distinctly different between on-and off-swirl regions, but both are consistent with compositionally equivalent global maria, labeled "On-swirl-like" and "Off-swirllike" (i.e., same FeO, OMAT, and TiO 2 ).

Figure 23 .
Figure 23.Map of Diviner regolith nighttime temperature anomaly (ΔT) at Chaplygin cold spot(Powell et al. 2023).We define our cold spot sample as the ∼167,000 LOLA data points inside the inner yellow circle for which ΔT is below the global 1st percentile (−3.65°C, assigned a color of white).The control sample is the ∼222,000 LOLA data points between the yellow circles.

Figure 24 .
Figure 24.The phase function of the Chaplygin cold spot is statistically indistinguishable from that of a surrounding control sample.The shaded regions show the interquartile range.

Figure 25 .
Figure 25.Phase functions for five silicic anomalies are similar to compositionally equivalent global highlands terrain.

Table 2
Coefficients for Piecewise Linear Fits of Phase Function Parameters in ±3°Equatorial Strip Using Fitting

Table 3
Coefficients from PLS Regression for Standardized Variables (To Show Relative Contributions) and Their Values in the Original Units for Each Variable

Table 4
Parameters for Phase Function Fits at Reiner Gamma Swirl

Table A1
Groupings of Tracks Used for the Dark Current and RADF/Counts Fits.