Spectrophotometric Modeling and Mapping of (101955) Bennu

Using hyperspectral data collected by OVIRS, the visible and infrared spectrometer on board the Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer (OSIRIS-REx) spacecraft, we modeled the global average spectrophotometric properties of the carbonaceous asteroid (101955) Bennu and mapped their variations. We restricted our analysis to 0.4–2.5 μm to avoid the wavelengths where thermal emission from the asteroid dominates (>2.5 μm). Bennu has global photometric properties typical of dark asteroids; we found a geometric albedo of 0.046 ± 0.007 and a linear phase slope of 0.024 ± 0.007 mag deg−1 at 0.55 μm. The average spectral slope of Bennu’s normal albedo is −0.0030 μm−1, and the phase-reddening parameter is 4.3 × 10−4 μm−1 deg−1, both over the spectral range of 0.5–2.0 μm. We produced normal albedo maps and phase slope maps at all spectral channels, from which we derived spectral slope and phase-reddening maps. Correlation analysis suggests that phase slope variations on Bennu are likely due to photometric roughness variation. A correlation between photometric and thermal roughness is evident, implying that the roughness of Bennu is self-similar on scales from tens of microns to meters. Our analysis reveals latitudinal trends in the spectral color slope and phase reddening on Bennu. The equatorial region appears to be redder than the global average, and the spectral slope decreases toward higher latitudes. Phase reddening on Bennu is relatively weak in the equatorial region and shows an asymmetry between the northern and southern hemispheres. We attributed the latitudinal trend to the geophysical conditions on Bennu that result in a global pattern of mass flow toward the equator.


Introduction
Near-Earth asteroid (101955) Bennu was the target of NASA's Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer (OSIRIS-REx) mission, which is returning a sample of carbonaceous regolith from this asteroid to Earth to understand the origin and evolution of volatiles and organics in the inner solar system (Lauretta et al. 2017. The spectral properties of Bennu measured prior to the spacecraft's arrival were confirmed by the OSIRIS-REx observations (Lauretta et al. 2019a). In particular, the data acquired using the OSIRIS-REx Visible and InfraRed Spectrometer (OVIRS; Reuter et al. 2018) in the visible to nearinfrared (VIS-NIR; 0.4-3.8 μm) spectral range confirmed that Bennu has a globally blue-sloped spectrum consistent with a B-type classification (Clark et al. 2011;Hergenrother et al. 2013;Simon et al. 2020b) and contains hydrated phyllosilicates (Hamilton et al. 2019), carbonates Simon et al. 2020b), and organics (Simon et al. 2020b). In addition, the OSIRIS-REx Thermal Emission Spectrometer (OTES; Christensen et al. 2018) in the mid-infrared (MIR) wavelengths detected the presence of magnetite (Hamilton et al. 2019), which has been confirmed by other OSIRIS-REx instruments (Lauretta et al. 2019a;Simon et al. 2020a).
Like other primitive carbonaceous (C-complex) asteroids, Bennu has a globally dark surface. Ground-based observations yielded a disk-integrated phase function of Bennu in the visible wavelengths with a slope of 0.04 mag deg −1 (Hergenrother et al. 2013). The visible geometric albedo was found to be 0.045 ± 0.005 based on ground-based photometry combined with the size derived from thermal infrared (Emery et al. 2014), consistent with the value of 0.044 ± 0.002 from the diskintegrated data collected by the OSIRIS-REx Camera Suite (OCAMS; Rizk et al. 2018) during the spacecraft's approach to Bennu . Using empirical photometric models to fit ground-based photometric data in the visible wavelengths, Takir et al. (2015) derived a geometric albedo of 0.047 ± 0.02 and a Bond albedo of 0.015 ± 0.005.
A more detailed analysis of phase reddening was performed by Fornasier et al. (2020) based on the measured spectral slopes at various spectral ranges and phase angles with the OVIRS spectra normalized at 0.55 μm. Their results showed a linear phase reddening with a coefficient of 4.4 × 10 −4 μm −1 deg −1 in the 0.55-2.5 μm range, consistent with Zou et al. (2021). Moreover, Hasselmann et al. (2021) studied the roughness of Bennu with a seminumerical statistical model aided by a raytracing technique using the OCAMS images and found an average rms slope of 27°with two possible components and possible specular reflection that accounts for about 2.6% of the total scattered light.
The disk-resolved data also revealed spatial variations in the spectral and photometric properties on Bennu. Although the average surface of Bennu is blue-sloped and dark, the asteroid is covered by rocks of diverse albedo and color, composition, and morphology but very well mixed at scales down to tens of centimeters (DellaGiustina et al. 2019(DellaGiustina et al. , 2021Lauretta et al. 2019a;Walsh et al. 2019;Kaplan et al. 2020;Simon et al. 2020b;Golish et al. 2021aGolish et al. , 2021b. Spectral clustering analysis with the OVIRS data indicated that the equatorial region on Bennu has a redder spectral slope than the global average (Barruci et al. 2020), consistent with the interpretation of ground-based spectral data in the 0.7-2.3 μm range (Binzel et al. 2015); however, this latitudinal pattern is not evident in the visible wavelengths . The spatial distributions of the spectral slope in the 0.55-2.5 and 0.45-0.55 μm ranges appear to be anticorrelated Fornasier et al. 2020). In addition, phase reddening on Bennu varies among boulders and craters, although no segregation is evident between these two geologic types, and no correlation with albedo, phase slope, or spectral slope is identified (Fornasier et al. 2020).
Part of this diversity is interpreted as primordial heterogeneity among Bennu's regolith particles Rozitis et al. 2020b). But it could also be influenced by spaceweathering processes such as those described in Pieters & Noble (2016), for example, solar wind and cosmic-ray particle bombardment, micrometeorite impact, UV radiation photolytic, solar heating, and others associated with exposure to the space environment. DellaGiustina et al. (2020), after examining the albedo and color of rocks and craters and their associated geological conditions, proposed that space weathering on Bennu causes a quick and strong bluing followed by a slow and weak reddening, ultimately resulting in Bennu's gently blue average spectral slope. Other processes, such as thermal fatigue (Lauretta et al. 2019b;Molaro et al. 2020aMolaro et al. , 2020b, meteoroid bombardment (Lauretta et al. 2019b;Ballouz et al. 2020;Bottke et al. 2020), subsurface water ice sublimation (Rozitis et al. 2020a), and dehydration (Lauretta et al. 2019b;Praet et al. 2021), could also alter the spectrophotometric properties of regolith particles.
OSIRIS-REx observations show a dynamic surface at various timescales, as indicated by the coexistence of young and old craters and characteristics of crater size frequency distribution (Bierhaus et al. 2019;Walsh et al. 2019;DellaGiustina et al. 2020). The spin rate of Bennu is accelerating owing to the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect Nolan et al. 2019), leading to a top-like shape  through ongoing mass movement (Jawin et al. 2020) and possibly other dynamic processes. The surface geopotential lows on Bennu are in the equatorial region , and material generally moves toward the equator (Jawin et al. 2020), potentially resulting in relatively recently exposed surfaces at mid-to-high-latitude regions. Particles are launched from Bennu's surface to space throughout its orbit (Lauretta et al. 2019b;Hergenrother et al. 2020), some of which fall back to the surface . The intersection of the rotational Roche lobe (where particles are dynamically trapped) with the surface of Bennu at about ±20°l atitude ) indicates a preferential collection of launched particles in the equatorial region (McMahon et al. 2020), which could be linked to the underdense equatorial bulge  and might contribute to relatively old crater-retention ages for the equatorial region. Large rocks distributed in the southern hemisphere (Jawin et al. 2020) appear to inhibit surface mass flow, leading to hemispherical differences in the distribution of slopes and roughness .
Taken together, these previous studies indicate that although Bennu has overall spectrophotometric properties typical of dark, C-complex asteroids, variations in these properties at the global and local scales are potentially related to distinct physical and geological processes. In this paper, we derive spectrophotometric properties and map their variations across the surface of Bennu by applying spectrophotometric modeling and mapping techniques to VIS-NIR spectral data collected by OVIRS, and we interpret the patterns in the context of geological and geophysical processes on Bennu's surface. Section 2 describes the data and the corresponding preparation for our analysis. Section 3 introduces the modeling and mapping process and results. In Section 4, we discuss the interpretations and implications of our results and present a correlation analysis between photometric properties and the surface geological and geophysical conditions on Bennu. Section 5 summarizes the main findings from our analysis.

Data Description
In our analysis, we used the same OVIRS data set (Reuter et al. 2019) as used by Zou et al. (2021) for their global photometric analysis. These data were collected during the Preliminary Survey (PS) and Detailed Survey (DS) phases of the OSIRIS-REx mission (see Lauretta et al. 2017 andLauretta et al. 2021 for an overview of mission sequences). The calculation of scattering angles (incidence angle, i; emission angle, e; and phase angle, α) used a facet-based shape model of Bennu derived from OCAMS images via stereophotoclinometry (SPC) with triangular facets of ∼0.8 m side length (SPC model v20; Barnouin et al. 2019Barnouin et al. , 2020. The shape model was first degraded to ∼3 m triangular facets by combining neighboring facets. Then, the scattering angles of each spectrum were calculated by averaging the angles corresponding to 100 points uniformly distributed within the OVIRS footprint. Given the much smaller facet size than the footprint size of the data we used (14-44 m), the averaging minimized the effects of the uncertainties in the shape model in the calculated scattering angles. We tested diskintegrated photometric modeling with the scattering angles calculated using different degraded resolutions of the shape model and found that the model results remain almost unaffected.
To avoid data close to the limb and terminator, where the uncertainties in scattering angles are relatively high and the performance of photometric models is usually poor, we discarded all of the spectra with incidence angle i > 70°or emission angle e > 70°in modeling. A total of 84,431 spectra were used in our analysis. Figure 1 shows the characteristics of the reflectance data. See Zou et al. (2021) for details about the data set we used.

Calibration and Scattering Geometry
The data were calibrated by the standard OVIRS calibration pipeline as described in Simon et al. (2018Simon et al. ( , 2021. The OVIRS instrument is a point spectrometer with five linear variable filters that cover the spectral range of 0.4-4.3 μm . The standard calibration of OVIRS data includes background subtraction, out-of-band (red leak) correction, radiometric calibration, and resampling to a standard wavelength grid with a spectral resolution of 0.2 μm from 0.4 to 2.4 μm and 0.5 μm from 2.4 to 4.3 μm . The resulting spectra are in units of intensity (W m −2 μm −1 sr −1 ). Finally, the thermal emission that dominates the spectrum from about 2.5 to 3.8 μm is removed (Simon et al. 2020b), and the reflectance spectrum is converted to a radiance factor (or I/F, where I is the scattered intensity and πF is incident solar flux) using a solar spectrum.
The absolute radiometric calibration is known to ±2% . However, two steps affect the relative (wavelength-to-wavelength) uncertainty in the calibrated reflectance spectra: thermal tail removal and out-of-band flux correction. Thermal tail removal is described in detail in Simon et al. (2020b). The processing is performed on a per-spectrum basis. After subtracting a solar spectrum that is reddened by a linear slope and scaled to the radiance at ∼2.1 μm, a singletemperature blackbody spectrum is fitted to the data to derive a temperature and scaling factor that match the thermal component of the data at 3.75 μm. Then, the best-fit blackbody spectrum is subtracted from the original spectrum. However, the best-fit scaling factors range between about 0.4 and 0.8, suggesting that the thermal tail is dominated by the hottest area inside the field of view (FOV) as characterized by the best-fit temperature, and the scaling factor approximately represents the fraction of the FOV that is covered by the hottest area. Our analysis suggested that the uncertainty of thermal removal is negligible within 2.0 μm and up to 2% within 2.5 μm, although it could be higher at longer wavelengths. We therefore focused our analysis at wavelengths <2.5 μm, and the results should be robust against the thermal removal uncertainty.
The correction of out-of-band flux involves an estimate of the flux leaked from long wavelengths to short-wavelength filters . It is removed by calculating the radiance at long wavelengths and applying correction coefficients (derived from ground test data) to the shorter wavelengths. The long-wavelength radiance is affected by surface temperature, detector sensitivity (which varies by detector temperature), and thermal contrast in the scene. The largest uncertainties are usually associated with spectra that have low overall intensity levels and contain high-contrast scenes, such as boundaries of shadows or different surface units.
In addition, changing the scene contrast as the FOV is scanned causes a more common artifact: discontinuities at the filter boundaries of 0.66, 1.1, 1.8, and 2.9 μm, etc. The spectra that display the most obvious segment boundaries are all from high-latitude areas and account for ∼1% of all spectra we used in our modeling. Although we discarded those spectra in our modeling, segment boundary discontinuities and all other factors varyingly contribute to both the absolute uncertainties and spectral shape, potentially with some systematic dependence on local solar time and latitude of each particular spectrum. Those uncertainties collectively affect our photometric modeling and show up as discontinuities in the fitted model parameter spectra.

Photometric Data Reduction
The data we used cover the whole surface of Bennu in almost all phase-angle stations, although strong correlations are evident between the latitude and the incidence and emission angles (Figures 1 and 2). The correlation between incidence angle and latitude is a result of the nearly 180°obliquity of Bennu. The correlation between emission angle and latitude reflects the geometry of data acquisition. The data we used were predominantly collected during the DS Equatorial Stations subphase, when the spacecraft stayed at several stations in the equatorial plane of Bennu that have fixed spacecraft-Bennu-Sun angles. The spacecraft pointed at the nadir direction of the target and scanned in a north-south direction to collect data, with Bennu rotating underneath to provide whole-surface coverage of the target. Therefore, the emission angle and latitude are correlated. This data acquisition geometry also resulted in discrete samplings in phase angle. During the DS Baseball Diamond subphase, the spacecraft was above and below Bennu's equatorial plane, providing some data points that broke the correlation between emission angle and latitude to some extent. The complicated topography on Bennu also helps expand the geometric coverages of our data in i and e.
To map the photometric parameters, we divided the surface of Bennu into small units and fitted photometric models to each unit independently using only the data within the corresponding unit. We adopted two division schemes for the mapping: a longitude-latitude grid with a 5°resolution in both directions (∼20 m at the equator) and a facet-based scheme using the SPC model v20  degraded to ∼12 m triangular facets. The reason for downgrading the resolution of the shape model is to match the footprint size of the OVIRS data that we used and keep the number of facets computationally reasonable in model fitting. The former simplifies the map projection of parameters and visual comparisons with other similarly projected maps but suffers from nonuniform areas in each grid cell and small numbers of data points in the cells at high latitude. The latter delivers nearly equal areas for all surface units, which facilitates the statistical analysis of parameter maps, as well as their 3D visualization wrapping onto the shape model of Bennu (Ferrone et al. 2021), but the resulting maps are hard to reproject in 2D. Therefore, we will use the longitude-latitude grid scheme to identify features in the maps and discuss their robustness (Section 3), and we use the facet-based scheme for 3D visualization and correlation analysis and interpretations (Section 4). The results from both schemes are compared for a consistency check.
The observation properties (number of spectra, illumination conditions, etc.) of the data for the longitude-latitude grid scheme are shown in Figure 2. To interpret spatial features in the maps of photometric parameters later, we need to compare them with these maps of observation properties to determine whether they might be modeling artifacts caused by the distributions of data points or scattering geometry. There is a slight hemispherical asymmetry in the number of spectra (Figure 2(a)), where slightly less spectra are available for the southern hemisphere. Overall, the minimum incidence angle distribution is symmetric with respect to the equator (Figure 2(b)). There is an incomplete coverage in emission angle south of −45°to −50°latitude, where no spectra with e < ∼50°are available (Figure 2(d)). The distributions of maximum incidence angle (Figure 2(c)) and emission angles (Figure 2(e)) are uniform over the whole surface with no systematic patterns visible. The phase-angle coverage (Figures 2(f) and (g)) shows that the region within ±40°l atitude receives the most complete coverage from <7°to >90°p hase, except for a small number of cells where the maximum phase angles are <∼40°. North of +40°latitude, the maximum phase angles are about 80°, whereas south of −40°latitude, the maximum phase angles are <∼40°. Within the ±40°latitude zone, the minimum phase angle (Figure 2(f)) shows some vertical bar-like patterns in the 6°-8°range. In two regions at high northern and southern latitudes, the minimum phase angle rises to about 10°. These patterns in phase-angle maps are all traceable to the observation geometry during various data collection campaigns. Finally, we rejected the spectra with extreme ranges of phase angles of <7°and >100°. All of the above will need to be considered when interpreting the derived photometric parameter maps.
The overall patterns in the data property maps for the facetbased division scheme are similar to what Figure 2 shows. However, the histograms of the number of data points in the units are different for the two schemes. For the longitudelatitude grid scheme, about 50% of cells contain >20 spectra, and about 20% of cells contain >40 spectra. In contrast, the facet-based scheme resulted in a more uniform distribution of data points over cells, with about 50% of cells containing >40 spectra, and about 20% of cells having >60 spectra.

Models and Fitting
We adopted the same models as in Li et al. (2019): the Akimov empirical model and the Hapke model. We briefly describe the models and the fitting approach in this section. Additional details can be found in Li et al. (2019).
For the empirical approach, we used the parameterless Akimov model (Shkuratov et al. 1999) where α, b, and l are the phase angle, photometric latitude, and photometric longitude, respectively. Then we fitted the phase function using an exponential model with two parameters, where A N is the normal albedo and β is the phase slope in units of mag deg −1 that characterizes the decrease of reflectance with respect to phase angle. The reflectance factor (RADF) of the surface is This model is easy to fit and has been demonstrated to be able to describe the scattering behaviors of many asteroids (e.g., Schröder et al. 2013Schröder et al. , 2017Longobardo et al. 2014Longobardo et al. , 2019Combe et al. 2015;Hasselmann et al. 2016;Li et al. 2019 (Hapke 1981(Hapke , 1984(Hapke , 1986) that we adopted here contains five parameters: single-scattering albedo (SSA), w; asymmetry factor, ζ, of the single-term Henyey-Greenstein (HG) function; photometric roughness, q; and the amplitude, B 0 , and width, h, of the shadow hiding opposition effect (SHOE). The functional form of the Hapke model that we adopted is where μ 0 and m 0 e are the cosines of local i and e corrected for roughness, respectively, and B SH is the SHOE. The form of B SH adopted here is the same as previously used in Li et al. (2004Li et al. ( , 2006; H (μ, w) is the Chandrasekhar H-function, for which we adopted the approximated form in Hapke (2002); S(q; i, e, α) is the correction for surface roughness that follows Hapke (1984); and p(ζ; α) is the single-particle phase function. With most of our data acquired at phase angles <95°, it is unnecessary to invoke a two-or three-parameter HG function to fit the single-particle phase function. We did not include anisotropic multiple scattering (Hapke 2002) in our modeling because multiple scattering is expected to be weak for the dark surface of Bennu. We did not include the coherent-backscattering opposition effect (CBOE) as introduced by Hapke (2002) because there are not sufficient data available at low phase angles to allow us to separate the CBOE from the SHOE and constrain their parameters, and also because the multiple scattering that causes CBOE is expected to be weak owing to Bennu's low albedo. We did not consider a porosity parameter (Hapke 2008) because it is impossible to separate the effect of porosity for a dark surface without sufficient data to constrain the opposition effect (Helfenstein & Shepard 2011). The geometric and Bond albedos follow their commonly adopted definitions and equations as given in Hapke (1981). Moreover, with a minimum phase angle of about 6°for our data, the SHOE parameters B 0 and h could not be reliably retrieved. Although the OSIRIS-REx spacecraft observed Bennu at phase angles down to nearly 0°during the Approach phase of the mission (Lauretta et al. 2017, the apparent size of the asteroid was smaller than the FOV of OVIRS, rendering the calibration of data unreliable for quantitative analysis. We had to use the disk-integrated phase function derived from the OCAMS MapCam data during Approach  to fit the Hapke model that we adopted here, and we derived B 0 = 2.06 and h = 0.11 as averaged over the four color bands (b, v, w, x) and the panchromatic filter. We then assumed that these two parameters are independent of wavelength over the OVIRS spectral range and held them fixed in fitting the global Hapke parameters. The values of B 0 and h affect the absolute uncertainty of the phase slope and SSA but not their band-to-band uncertainty. The effect on the roughness parameter is also insignificant because that parameter is mostly determined by the limb-darkening characteristics and the shading effect from local topography.
The model fitting follows a least-squares approach, which minimizes the squared sum S of the model fitting residual, where r i is the measured RADF, r i,model is the modeled RADF, and the sum is over all n data points. Here S is minimized to derive the best-fit parameters using the Levenberg-Marquardt algorithm with a constrained search space for the model parameters (Moré 1978, Markwardt 2009). We also define the relative rms (rel. rms) as where S is the squared sum in Equation (6), and r is the average measured RADF of all data points. The rel. rms is normalized by the number of data points and average measurements and allows for comparison of the model fits to different data sets to some extent. The OVIRS spectra have 1393 spectral channels, and we fitted all channels independently. For photometric parameter mapping, we independently fitted all of the grid cells that contain more than 10 spectra.
We followed the approach of Li et al. (2013cLi et al. ( , 2019 to estimate the uncertainties of the best-fit parameters, taking into account the possible correlation between model parameters. In particular, to estimate the uncertainty of a particular parameter, we fit the model with that parameter fixed in a series of values within a range surrounding the best-fit value to derive the corresponding minimum S values. We define the 1σ uncertainty of this parameter as the range where S is less than 2× the bestfit S. The uncertainties estimated from this approach should be considered systematic error because they take into account all sources of uncertainty, such as the measurement error (both absolute and statistical), the uncertainties in the scattering angles caused by shape model uncertainties, the imperfection or incompleteness of the model, and the correlation between model parameters. Therefore, they are likely the most conservative uncertainty estimate that should be used when interpreting the photometric parameters in terms of the absolute physical properties of the surface and when comparing with other objects.
On the other hand, when analyzing the spectral features of model parameters and the spatial patterns of parameter maps, what matters is the relative channel-to-channel and location-tolocation uncertainties, respectively. The relative uncertainties are dominated by the statistical uncertainties of the data but minimally affected by factors such as the absolute calibration error or the imperfection of model. The correlation between model parameters also affects the relative uncertainty to some extent depending on the characteristics of the data and model. A good measure of the relative uncertainties is the noise in the parameter spectra and maps that can be quantified by rolling standard deviations within a certain width in the spectra and maps. We will discuss the robustness of the spectral and map features primarily based on the relative noise and take the systematic or absolute uncertainties of the parameters as a reference.

Global Photometric Model
Before deriving the maps of the photometric parameters, we fitted both models to all spectra to derive the global average model parameters. The global modeling serves as a baseline for a consistency check of the mapping and a basis for the analysis and interpretations of maps.

Akimov Empirical Model
The goodness of fit for the Akimov empirical model fitting at 0.55 μm is shown in Figure 3. Figure 4 shows the best-fit global parameter spectra. Overall, the rel. rms for global modeling is about 6%, which is comparable with previous modeling results of other asteroids and cometary nuclei based on multispectral or hyperspectral data collected during spacecraft visits (e.g., Li et al. 2009Li et al. , 2013aLi et al. , 2013bLi et al. , 2013cLi et al. , 2019Schröder et al. 2013Schröder et al. , 2017Ciarniello et al. 2017), suggesting satisfactory modeling results. No systematic trend with respect to any scattering angle is evident for the best-fit model (Figures 3(b)-(d)), implying that our fitted model is able to satisfactorily describe the photometric behavior of Bennu's surface.
In terms of wavelength trend, the rel. rms remains stable at wavelengths shorter than about 2.5 μm, then starts to increase dramatically and forms a bump in the range of 2.7-3.5 μm (Figure 4(a)), suggesting a poorer fit in this spectral range. This degradation of fit is likely due to the potential imperfection of thermal tail removal (Section 2.2). There are also some small discontinuities at about 0.65, 1.1, 1.8, and 2.3 μm that are likely associated with segment boundaries in the OVIRS data, suggesting varying calibration qualities over the segments and sudden changes at the boundaries. The noise in the parameter spectra appears to be very low compared to the overall spectral shapes and features, suggesting small band-toband uncertainty and reliable spectral shapes.
The fitted normal albedo spectrum (Figure 4(b)) shows an overall blue-sloped shape and a 2.7 μm hydroxyl absorption feature. The normal albedo at 0.55 μm is 0.046 with an error of about 14% (see Section 3.2.3 for error analysis). Our empirical phase function model does not include the opposition effect, which is about 10% . Therefore, the true normal albedo derived from our data, taking into account the opposition effect, would be about 0.051. For a surface that scatters light following the parameterless Akimov disk-function model, normal and geometric albedo are equal. Because all previous ground-based analyses (Hergenrother et al. 2013;Emery et al. 2014;Takir et al. 2015) and the early measurements from OSIRIS-REx OCAMS images  took into account the opposition effect, our modeling resulted in about 10% higher geometric albedo than those values, although it is still consistent within their measurement uncertainties. This apparent discrepancy is in fact consistent with the absolutely calibrated flux level of OVIRS being ∼8%-10% higher than those of OCAMS and OTES (Rozitis et al. 2020b;Zou et al. 2021).
The phase slope at 0.55 μm is fitted to be 0.024 mag deg −1 . This value translates to 0.033 mag deg −1 for the disk-integrated phase function by accounting for the decreasing visible and illuminated fraction with respect to phase angle, and it is shallower than the value of 0.039 mag deg −1 derived from both ground-based data and OCAMS imaging data acquired during Approach (Hergenrother et al. 2013Takir et al. 2015). The linear spectral slope between 0.5 and 2.0 μm is fitted to be −0.0030 μm −1 . Fornasier et al. (2020) reported a linear spectral slope of −0.043 μm −1 between 0.55 and 2.5 μm based on the global spectrum of Bennu at about 8°phase angle normalized at 0.55 μm. If we use the same wavelength range as used by Fornasier et al. (2020) and normalize the fitted normal albedo spectrum at 0.55 μm, the linear spectral slope is −0.047 μm −1 . Further, using the phase-reddening parameters of 4.3 × 10 −4 (μm deg) −1 fitted from our modeling, we derive a slope of −0.044 μm −1 at 8°phase angle, consistent with the Fornasier et al. (2020) value. The spectral shape of the longwavelength side of the 2.7 μm feature appears to be slightly different from OVIRS spectra (Hamilton et al. 2019). Because of the increased model rms beyond 2.7 μm, as also evidenced by the increased noise in the normal albedo spectrum (Figure 4), the band shape in the normal albedo spectrum may not be reliable.
The phase slope parameter shows an overall decrease with wavelength (Figure 4(c)), consistent with phase reddening. The phase-reddening coefficient derived as the spectral slope of the phase slope parameter is 4.3 × 10 −4 (μm deg) −1 , consistent with the Zou et al. (2021) result using a similar approach and the Fornasier et al. (2020) result using a different approach. A sudden drop in phase slope appears on the short-wavelength side of the 2.7 μm feature and forms a steplike shape inside the band (Figure 4), indicating decreasing band depth with respect to phase angle and minimal change in the shape of the band on the long-wavelength side, consistent with the observations of Fornasier et al. (2020). The segment boundaries at 1.1 and 1.8 μm are clearly visible in the phase slope spectrum, although not apparent in the normal albedo spectrum. This is evidence that the uncertainties in the OVIRS data calibration, although not visually significant in most of the spectra, propagate through the photometric modeling and show themselves in the model parameter spectra. Finally, the increased noise in the phase slope parameter spectrum, another indication of degraded model fit quality in this spectral range, suggests that the spectral shape here is unreliable and should not be interpreted.

Hapke Model
The goodness-of-fit plot for our global Hapke model at 0.55 μm is shown in Figure 5, and the best-fit parameter spectra are shown in Figure 6. Similar to the goodness of fit of the Akimov empirical modeling, there is no systematic trend evident (Figures 5(b)-(d)), suggesting that the Hapke model we adopted was able to describe the photometric behavior of Bennu. The rel. rms is slightly lower than the Akimov empirical model fitting, indicating a slightly better model fit to the data. All other goodness-of-fit features are similar between the two models.
The parameter spectra for the Hapke model also show similar characteristics to the Akimov empirical model parameter spectra, where, shortward of the 2.7 μm absorption feature, the spectra are flat or roughly linear, and a feature or a sudden jump appears at about 2.7 μm with increased noise. Again, the model parameters at or beyond 2.7 μm may not be reliable for interpretations.
The characteristics of the three fitted Hapke parameters (w, ζ, and q) and the two derived parameters (geometric albedo, p v , and Bond albedo, A) are as follows.
1. The SSA spectrum is mostly flat at <2.5 μm, shows the 2.7 μm hydroxyl absorption feature, and then becomes very noisy beyond 3.2 μm (Figure 6(b)). The SSA is about 0.055 before the 2.7 μm feature and 0.047 near the band minimum. 2. The asymmetry factor ζ shows a linear spectrum before the 2.7 μm feature, increasing from about −0.27 at 0.5 μm to −0.25 at 2.5 μm (Figure 6(c)). A small bump centered near the 2.7 μm absorption center has a maximum ζ of −0.22. Outside of the band, the ζ spectrum appears to continue the linear trend. Because less negative values mean less backward scattering, given the nearly flat q spectrum (see next bullet item), the increasing trend of ζ is consistent with phase reddening on Bennu, and the 2.7 μm bump is consistent with weaker hydroxyl absorptions at higher phase angles. 3. The roughness parameter q remains in a narrow range between 22°and 25°(Figure 6(d)), which is consistent with our expectation that, as a geometric parameter, it is not expected to depend on wavelength. We do not attempt to interpret the spectral shape of q. 4. The spectra of the geometric (Figure 6 In contrast, the Bond albedo depends on the overall shape of the phase function, and the values are therefore more reliably determined than the geometric albedo. Similar to all other parameter spectra, the values are not reliable beyond about 2.5 μm. 5. At 0.55 μm, the modeled geometric albedo of Bennu is about 0.053, and the Bond albedo is about 0.021. This value is consistent with the value derived from the Akimov empirical model when considering the opposition effect but higher than those previously derived from ground-based observations and the OCAMS Approach data , again due to the difference in the absolute flux calibration of OVIRS and OCAMS (Rozitis et al. 2020b;Zou et al. 2021).

Global Model Uncertainty Analysis
The uncertainty analysis of the global model parameters follows the general procedure as described in Li et al. (2013cLi et al. ( , 2019 and outlined in Section 3.1. We considered the uncertainties of the best-fit parameters at three wavelengths-0.55, 3.0, and 3.7 μm-where the rel. rms indicates the best fitting quality, the worst fitting quality, and where the best-fit parameter spectra are the noisiest, respectively. The uncertainty analysis for the Akimov empirical model parameters is shown in Figure 7. The 1σ errors for the best-fit normal albedo and phase slope are ∼7.8% and ∼19%, Figure 4. Global averaged model parameters and rms for the Akimov model. The model rms dramatically increases from about 2.7 μm and could be associated with the increased uncertainty in reflectance data in this spectral range due to thermal tail removal. The noise in the rms and parameter spectra becomes excessively high beyond about 3.3 μm, owing to the low signal in reflectance after thermal tail removal. respectively, at 0.55 μm (Figures 7(a) and (d)). For the worst fitting quality at about 3.0 μm, the relative errors are about 13% and 33%, respectively (Figures 7(b) and (e)). At 3.7 μm, although the fitted parameter spectra are noisy, the uncertainties seem to be determined by the rel. rms and are just slightly higher than those at 0.55 μm (Figures 7(c) and (f)). The lower bounds of both normal albedo and phase slope are slightly better constrained than the higher bounds. As discussed earlier, this uncertainty estimate should be considered systematic in that it includes all factors that affect the modeling, although it does not include the absolute radiometric calibration uncertainty of the OVIRS instrument of about 2% ). On the other hand, the noise in the spectra themselves offers a more direct way to indicate the band-to-band uncertainty. The increased noise in the parameter spectra at beyond 3.0 μm is less than the uncertainty derived above. Therefore, the wavelength dependence of the best-fit parameters (Figure 4) should still be considered valid, despite the fact that the ranges of variations for both spectra are compatible with or sometimes higher than the model uncertainties.
In contrast, the error analysis for the fitted Hapke parameters (Figure 8) suggested that those parameters were not well constrained in an absolute sense. The lower bounds of SSA and ζ were about 0.044 and −0.35, respectively, for the best case at 0.55 μm (Figures 8(a) and (d)). However, their higher bounds were poorly constrained to about 0.09 and −0.14, respectively. The situation is worse for the other two bands at 3.0 and 3.7 μm (Figures 8(b) and (e)). The lower bound of the roughness parameter is not constrained at all by the data we used in any bands, while the higher bounds are unconstrained for the worst case and between 45°and 50°for other cases (Figures 8(g), (h), and (j)). Any value of roughness between 0°and 50°could fit the data without increasing the best-fit S in a statistically significant way, meaning that no values of q in the possible range can be rejected by the fit.
The reason for the poor constraint on the Hapke parameters is likely due to the distinctive surface texture of Bennu. The OCAMS images of Bennu reveal a rugged surface almost fully covered by rocks of sizes comparable to or larger than the image pixel scale that have highly irregular and angular shapes and flat facets on the surface (e.g., Lauretta et al. 2019aLauretta et al. , 2019b. There are no extensive smooth areas on Bennu that appear to be uniformly covered by fine-particulate materials like the ponded areas on Eros (Cheng et al. 2002) and Itokawa (Fujiwara et al. 2006). Bennu also appears to be distinct from the Hayabusa2 mission's target, asteroid (162173) Ryugu, which, although also almost entirely devoid of smooth regions, has smaller and much less angular rocks than those on Bennu These factors combined cause a large scatter for Bennu's surface when fitting to the commonly adopted light-scattering models. Therefore, the Hapke roughness parameter, which dominates the reflectance distribution with respect to local topography, cannot fully account for the reflectance variations and is poorly constrained. As a consequence, the asymmetry factor, which is primarily fitted from the phase slope, cannot be well constrained because the roughness parameter also affects the phase slope. Ultimately, the strong correlation between SSA and asymmetry factor leads to poor constraints on SSA.
The large uncertainties of the modeled Hapke parameters hinder their physical interpretations. Therefore, we made no attempt to interpret these Hapke parameters and compare with other objects. On the other hand, these uncertainties are systematic, and the bandto-band uncertainty of the fitted parameter spectra is best estimated from the noise in the spectra (Figure 6). Those noises appear to be much smaller than both the systematic uncertainties we derived here and the spectral features, except for spectral bands >3.5 μm. Therefore, the analysis and interpretations of the spectral features in the fitted parameters should still be valid.

Akimov Empirical Model
By fitting the photometric models to each surface unit following both the longitude-latitude and facet-based division scheme, we derive map cubes of all photometric parameters and the rel. rms. Figure 9 shows the maps at 0.55 μm derived using the Akimov empirical model in both schemes. The rel. rms is between 1% and 6% with a random pattern, showing satisfactory goodness of fit for all cells (Figure 9(a)). Inspection of the goodness-of-fit plot (similar to Figure 3) for some randomly selected cells suggests no systematic trend in the model with respect to scattering geometry. As noted previously, when analyzing the spatial distribution of parameter values in the map, the point-to-point relative uncertainty should be considered, instead of the systematic uncertainty as described earlier for global models. The relative uncertainty in the maps can be represented by the noise and quantified by the rolling standard deviation map. Calculation shows that for the A 0 map at 0.55 μm (Figures 9(b) and (d)), the rolling standard deviation with a box size of three to five points on a side is typically <2% of the global average and <15% of the range of values in the map, except for a few small areas in the darkest regions that can be up to twice as high. For β at 0.55 μm (Figures 9(c) and (e)), those values are about 4% and 10%, respectively, and almost uniform over the map. Moreover, the spatial features in A 0 and β are overall consistent with only slight and continuous change across wavelength bands. Therefore, the features in the maps are reliable.
The normal albedo map that we derived (Figures 9(b) and (c)) is consistent with the albedo maps previously derived from OCAMS (Golish et al. 2021b) and OVIRS (Zou et al. 2021) data using the traditional photometric correction approach. The phase slope map (Figures 9(c) and (e)) shows some similar pattern in the dark areas in the southern hemisphere, such as the dark boulder Roc Saxum and the Tlanuwa Regio, where low albedo appears to correlate with higher (steeper) phase slopes. But no such correlation is evident in the northern hemisphere and for areas with average-to-high albedos. The lack of patterns in the rms map (Figure 9(a)) and the evident patterns in the parameter maps (Figures 9(b)-(e)) suggest that the distributions  Figure 4 for the Akimov model, the increase in model rms and noise in the rms and parameter spectra beyond about 2.7 μm are due to relatively high uncertainty and a low signal level after thermal tail removal. of normal albedo and phase slope are unlikely to be due to modeling artifacts. No patterns similar to those in the observation property maps (Figure 2) are seen either. Moreover, the maps derived using the facet-based scheme show similar patterns. Therefore, we have confidence that the albedo and phase slope patterns we observe are real.
To compare the mapping results with the global modeling results, we calculated the mean and median for each spectral channel over the whole map from both schemes and compared with the best-fit parameter spectra derived from global modeling ( Figure 10). The mean and median parameter spectra are within 0.5% of each other and 1% of the global modeling results for both division schemes, showing good agreement between global model and model parameter mapping.
Each derived map is a spectral cube with 1393 spectral channels along the spectral dimension, allowing us to study the spectral behavior of the photometric parameters. Based on the nearly linear spectral shapes of the normal albedo and phase slope parameters before the 2.7 μm absorption (Figures 4 and  10), we fitted the linear slopes S of the two parameters in various spectral bands between 0.44 and 2.3 μm, defined as where λ 1 and λ 2 are the wavelengths of the two boundaries and y(λ) is the fitted normal albedo or phase slope at wavelength λ.
The slope that we calculated does not normalize to a specific wavelength but rather has the unit of the corresponding quantity per wavelength. We stop at 2.0 μm to avoid the slight curvature in the spectra of both parameters starting from about this wavelength. The noise in the model parameter spectra is high beyond the 2.7 μm band due to thermal tail removal, as discussed in Section 2.2. We therefore do not include the 2.7 μm band and longer wavelengths in our analysis.
In addition, as suggested by DellaGiustina et al. (2020), the color properties of Bennu's surface depend on spectral range in the visible wavelengths. We divided the OVIRS spectral range into several segments and fitted the spectral slopes in the segments separately. As discussed in Sections 2.2 and 3.2, the filter segment boundaries are visible in the model parameter spectra. We aligned our spectral segments with the filter segments and discarded the channels within 0.02 μm from the segment boundaries at 0.66, 1.09, and 1.8 μm in order to minimize the effect of this calibration uncertainty. The boundary of the last spectral segment at 2.3 μm is determined by the apparent curvature in the spectra entering the 2.7 μm absorption. Some small spectral features of up to a few percent in strength centered near 1.05, 1.8, and 2.3 μm appear in the global normal albedo spectrum (Simon et al. 2020a), but they are weak compared to the overall spectral shape and should not appreciably affect the spectral slopes that we fitted across each spectral segment and the full range from 0.5 to 2.0 μm. Figure 11 shows the spectral slope maps derived from the normal albedo cube in various spectral segments. Figure 12 shows the spatial trend of the spectral slopes of Bennu's surface with respect to longitude and latitude. The spectral slope maps based on the facet-based scheme are consistent with Figures 11 and 12. Bennu's surface shows an overall blue-sloped spectrum, in agreement with previous results (Clark et al. 2011;DellaGiustina et al. 2020;Simon et al. 2020a). In the spectral range of 0.5-2.0 μm, the equatorial belt of Bennu appears to be distinctly redder than the rest of the surface (Figure 12(a)). Other relatively red areas include, for example, the OSIRIS-REx mission's primary sample collection site, Nightingale; the backup sample site, Osprey; and the dark Roc Saxum and Tlanuwa Regio. The spectral slope in the range 0.44-0.64 μm (Figure 11(b)) shows a general anticorrelation with the spectral slope of 0.5-2.0 μm (Figure 11(a)), as also noticed by Fornasier et al. (2020). The relatively blue equatorial region in the 0.44-0.64 μm range has a slightly different spatial distribution from the redder equatorial belt in the 0.5-2.0 μm range, indicating deviation from the trend for some areas. Moving to longer wavelengths, the 1.11-1.78 μm spectral slope map (Figure 11(d)) appears to dominate the spectral slope of 0.5-2.0 μm, likely because it is the widest range of the four segments. The map in the last spectral segment, 1.82-2.3 μm Figure 7. Error analysis of the best-fit photometric parameters for the Akimov empirical model. As labeled in each panel, the top row is for normal albedo, and the bottom row is for the phase slope parameter; the columns from left to right correspond to three bands at 0.55, 3.0, and 3.7 μm. The abscissae are parameter values, and the ordinates are the squared sum of the residual, S. The two horizontal dotted lines mark the minimum S and twice that value. The vertical dashed line marks the bestfit value for the corresponding parameters.
( Figure 11(e)), appears to be very noisy, with some patterns similar to those in the observation properties maps (Figure 2(f)), such as the horizontal band at +40°latitude and the vertical linear pattern at low latitude close to 270°longitude. The maps derived from this spectral segment should be excluded from our interpretation. Despite the apparent latitudinal trend of spectral colors, there is no evidence for a longitudinal trend (Figure 12(b)).
Interestingly, equatorial reddening is not observed in the visible wavelengths from the color ratio maps produced from OCAMS multiband images; instead, these maps suggest a spectral bluing in the equatorial region in the 0.55/0.47 μm ratio and no latitudinal trend in the 0.85/0.55 μm ratio (DellaGiustina et al. 2020). For a consistency check, we generated OVIRS color slope maps in those two spectral ranges ( Figure 13). Our maps appear to be consistent with the corresponding OCAMS color ratio maps, although the OVIRS maps are noisier and the spatial resolution is lower. No latitudinal trend is evident in the 0.55-0.85 μm spectral slope map. This agreement suggests that our photometric modeling and mapping approach did not introduce latitudinal color trends, and the spatial distribution of the spectral slope on Bennu is wavelength-dependent. Figure 14 shows the spectral slope maps of the phase slope parameter. The phase slope measures how quickly the reflectance decreases with respect to the phase angle. If the phase slope becomes smaller (shallower) with respect to wavelength, this means an increased spectral slope with respect to the phase angle, that is, phase reddening. Therefore, the spectral slope of the phase slope parameter quantifies phase reddening, with negative values corresponding to phase reddening and positive values to phase bluing. These maps are thus equivalently phase-reddening maps. Figure 15 shows the longitudinal and latitudinal trends of phase reddening. The surface of Bennu has an overall phase reddening, as reported by many other studies using OSIRIS-REx data (DellaGiustina et al. 2019;Hergenrother et al. 2019;Fornasier et al. 2020;Golish et al. 2021a;Zou et al. 2021). The reddening is strongest in the shortest-wavelength segment (Figure 14(b)) and decreases with respect to wavelength. Beyond 1.8 μm, phase reddening almost disappears or flips to slight bluing (Figure 14(e)). The spectral trend of phase reddening on Bennu from our analysis is consistent with the findings of Fornasier et al. (2020).
The spatial distribution of phase reddening on Bennu does not show obvious localized patterns above the noise, as also shown by Fornasier et al.'s (2020) analysis with selected boulders and craters. On the other hand, the northern hemisphere appears to have a slight systematically higher phase reddening than the southern hemisphere, and the phase reddening in the 0.5-2.0 μm spectral range might be weaker along the equatorial band and become more enhanced toward high latitudes (Figure 14(a)). No other latitudinal or longitudinal trend is visible. Some potential artifacts start to appear in the 1.11-1.78 μm range, such as the horizontal strips at ∼±45°latitude (Figure 14(d)). The single-line horizontal strip along the equator could be an artifact associated with the number of data point distributions (Figure 2(a)). In the 1.82-2.3 μm segment (Figure 14(e)), the northern hemisphere shows stronger phase bluing than the southern hemisphere. But again, the noise and artifacts in this spectral range may mean that this pattern is not reliable.

Hapke Model
The photometric model mapping results based on the Hapke model are shown in Figure 16. Similar to the global modeling, we fixed the opposition parameters B 0 = 2.06 and h = 0.11. The rel. rms map is similar to that derived from the Akimov empirical model (Figure 9) in distribution, and the absolute level is slightly lower, implying slightly better model fitting.
On the other hand, the maps of the three fitted Hapke parameters-the SSA (Figure 16(d)), asymmetry factor (Figure 16(b)), and roughness parameters (Figure 16(c))-are all noisy. The derived geometric albedo map is noisy as well, but the Bond albedo map appears to be much smoother. The noise in the Hapke model maps results from the relatively high model uncertainty for our Hapke model fitting, due primarily to the small number of data points in each grid cell compared to previous similar work (Li et al. 2019). The noise level naturally provides a rough estimate of the statistical uncertainties of the modeled parameters.
Despite the noise in the maps, some patterns are still visible. All three albedo maps show similar bright and dark patterns similar to the normal albedo map derived from the Akimov empirical model (Figures 9(b) and (c)). The reason for the much lower level of noise in the Bond albedo map (Figure 16(f)) than in the other two albedo maps (Figures 16(d) and (e)) is that the Bond albedo depends on the overall shape of the observed phase function and is therefore insensitive to the model fitting process, whereas both the SSA and the geometric albedo are highly model-dependent. The noise level in the SSA map is estimated to be about 0.01, about 25% of the respective averages. The asymmetry factor map has a similar pattern to the phase slope parameter map from the Akimov empirical model (Figure 9(d)), where more negative values, which imply stronger backscattering, correspond to steeper phase slopes. The roughness map (Figure 16(c)) is full of random noise, and it is hard to reliably discern any patterns. The standard deviation of the roughness map is about 9°, which is a reasonable estimate of the relative uncertainty for the fitted roughness parameter.
Similar to the Akimov empirical model results, we plotted the mean and median Hapke parameter spectra over the whole surface, together with the global model parameter spectra in Figure 17. The differences between the mean and median are much higher than those for the Akimov empirical model, again reflecting the high uncertainty, both systematic and statistic, of the Hapke modeling results. A small number of model parameters that are a few sigma away from the average would drag away the mean from the median. Therefore, it can be inferred that the S space of the Hapke model fitting is flat around the minimum, resulting in large uncertainties and unstable best-fit parameters across the surface, causing high noise in the maps. Again, for the Bond albedo, the good agreement between mean, median, and global model results is evidence that the Bond albedo is insensitive to the model fitting process and much more reliable than other parameters.
The spectral analysis based on the Hapke parameter maps resulted in very noisy maps for the spectral slope maps and phase-reddening maps, which thus are not of value for further interpretations. We therefore did not perform Hapke model fitting to the facet-based division scheme data, and we will base

Discussion
To decipher what the normal albedo and phase slope maps indicate about Bennu, we analyzed correlations within pairs of various maps. Included in our analysis are the normal albedo map, phase slope map, spectral color map, and phase-reddening map in the range 0.5-2.0 μm, as well as the thermal inertia and thermal roughness maps (Rozitis et al. 2020b). The correlation analysis is based on our facet-based division scheme using the shape model with a triangular facet side length of ∼12 m and extrapolated to a triangular facet side length of ∼6 m in order to match the resolution of previously generated thermal maps (Rozitis et al. 2020b). The use of the facet-based division scheme in the quantitative analysis ensures that the correlation analysis is not biased by map projections. The scatter plots relating normal albedo to other quantities and the corresponding correlation coefficients, R, are shown in Figure 18. In addition, to study the correspondence of spectrophotometric properties with geological features on Bennu, the maps are wrapped around the shape model of Bennu and overlaid on the OCAMS base map (Bennett et al. 2020), as shown in Figure 19. and phase slope (c) of Bennu based on the Akimov empirical model from the two division schemes. The gray lines are the parameters derived from global modeling (Section 3.2.1). The rel. rms spectrum from global modeling (not shown) is much higher than the rel. rms shown here. The blue and orange curves are based on the two schemes as noted in the legend. The solid lines are medians, and the dotted lines are means. For the longitude-latitude maps, the calculation of mean and median is based on reprojection to the equal-area sinusoidal projection. The yellow shaded area marks the spectral range from 0.5 to 2.0 μm; the gray shaded areas mark the spectral regions excluded from our analysis. The four spectral ranges that we used to fit spectral slopes are 0.44-0.64, 0.68-1.07, 1.11-1.78, and 1.82-2.3 μm. Figure 11. Normal reflectance slope maps of Bennu (in units of μm −1 ) in various spectral ranges based on the Akimov empirical model.

Surface Roughness
To interpret the albedo and phase slope maps, the first question we asked was, what dominates the phase slope map, surface roughness, or single-particle scattering characteristics? Photometric parameter mapping with the Hapke model could not provide a reliable answer to this question. But we can gain some insight from the correlation analysis between the phase slope and albedo maps. Figure 18(d) shows a non-or very weak correlation between albedo and phase slope for Bennu on a global scale, although some correlation could exist in local regions, such as those near the dark boulder Roc Saxum and the dark Tlanuwa Regio, as discussed earlier (Figures 9(b) and (c) and 19). A general correlation exists between albedo and phase Bennu's color slope shows a latitudinal trend, where the equatorial region appears to be relatively red at wavelengths longer than 0.6 μm but bluer in the wavelengths ranging from 0.44 to 0.64 μm. No such longitudinal dependence is evident in the full spectral range of 0.44-2.3 μm. Figure 13. Same as Figure 11, but in the spectral ranges that correspond to the OCAMS color ratio maps presented by DellaGiustina et al. (2020, their Figures  1(c) and (d)). Both maps are consistent with the corresponding OCAMS color ratio maps. Figure 14. Spectral slopes from phase slope maps of Bennu (in units of μm −1 deg −1 ) as a proxy for phase reddening, derived from the Akimov empirical model. slope in asteroids, where bright objects, or bright materials on a particular object, usually have shallower phase slopes than dark objects or materials, and vice versa (e.g., Belskaya & Shevchenko 2000;Schröder et al. 2013). The lack of a global correlation between albedo and phase slope on Bennu indicates that the phase slope map is likely not dominated by albedo but rather by photometric roughness. We can thus use the phase slope map as a proxy for photometric roughness. Both direct comparison ( Figure 20) and correlation analysis (Figure 18(f)) indicate that photometric and thermal roughness are highly correlated. The thermal roughness map is in good agreement with the rms slope roughness calculated from the global 20 cm resolution shape model (Rozitis et al. 2020b), which was found to be sensitive to roughness features as small as about 20 cm and as large as the resolution of the shape model used for the global thermal analysis (6 m). The thermophysical behavior (i.e., diurnal surface temperature distribution) is not strongly affected by roughness features smaller than 20 cm due to lateral heat conduction that minimizes any temperature deviations that would arise from the distribution of surface normals and from shadowing. The thermal roughness map qualitatively correlates with the distribution of boulder fields where the boulders are meter-scale and greater in size. The physical interpretation of photometric roughness, on the other hand, is not fully established but probably dominated by roughness at the smallest scale that still casts shadows (Shepard & Campbell 1998). In the case of our data, this scale is about 10-100 μm. The correlation between thermal and photometric roughness therefore implies that the roughness on Bennu is selfsimilar on scales from tens of microns to meters. OCAMS images of Bennu (e.g., DellaGiustina et al. 2019;Walsh et al. 2019) showed that the surface roughness at the centimeter scale and larger is dominated by rocks of centimeters to meters that are ubiquitous on its surface. If similar roughness exists at a scale of tens of microns, then it could be caused either by the existence of particles of that size, the roughness on the surfaces of rocks and large particles, or both.
Conflicting evidence exists for the presence of fineparticulate material on Bennu. On one hand, the low surface gravity of Bennu does not favor fine particles. Electrostatic levitation can preferentially loft fine-particulate dust from the surface of bodies such as Bennu (Hartzell & Scheeres 2013;Hartzell et al. 2019). Once ejected, it is harder for small particles to fall back to Bennu than for large particles (McMahon et al. 2020). The apparent flattening of the size frequency distribution of particles actively ejected from Bennu's surface (Hergenrother et al. 2020) could continue into micron-sized particles and might reflect a relative deficit of very small particles on the surface. On the other hand, on the basis of laboratory experiments, micron-sized particles could explain the color of Bennu's surface in comparison with laboratory samples (Sen et al. 2021). Spectral mixing models using laboratory-measured spectra of potential constituent minerals of Bennu's surface also require micron-sized grains in order to fit the absorption feature at 2.7 μm in Bennu's spectra (Merlin et al. 2021). Similarly, fine dust is also required to fit the OTES spectra for some regions on Bennu (Hamilton et al. 2020). Additionally, the results from thermal analysis are consistent with fine particles in localized regions (Rozitis et al. 2020b).
If the photometric roughness of Bennu is controlled by roughness at the tens of microns scale on the surfaces of rocks, then variations in the roughness are probably caused by heterogeneous compositions and/or physical properties within individual boulders and between different boulders (Della-Giustina et al. 2020, as well as the breakdown of rocks due to thermal fatigue (Molaro et al. 2020a(Molaro et al. , 2020a and impactinducted fracturing (Ballouz et al. 2020). Dark rocks on Bennu have rougher, more undulating surface textures than bright rocks . This is probably the reason for the weak correlation between high photometric roughness and low albedo in some regions in the southern hemisphere. If, on the other hand, the photometric roughness of Bennu is caused by particles of tens of microns in size, then highroughness regions would be correlated with a relatively high abundance of those particles. The weak correlation between high-roughness and dark regions is probably owing to the fact that dark boulders are mechanically weak and more prone to break, producing more small particles than their brighter counterparts. But if the apparent depletion of fine regolith on the surface of Bennu is real, then this also means that either the breakdown of particles stops at tens of microns or some processes must be continuously removing finer particles efficiently. Hasselmann et al. (2021), who studied Bennu's roughness by applying a statistical-numerical model (Van Ginneken et al. 1998) to the four-color OCAMS data, estimated an average roughness slope of  - +  20 . 5 . 1 . , very similar to our Hapke roughness slope of ∼23°at 1 μm. However, their Markov Chain Monte Carlo analysis also detected a second solution for the roughness slope of  - +  11 . 6 .
3 . , indicating that at subcentimeter size scales, the surface roughness slope distribution may be non-Gaussian, which is a general assumption of the Hapke shadowing function, or "populated" by two distinct roughness populations, one much less rough than the other. If, at the optical regime, a higher roughness slope corresponds to a smaller size scale through surface or particle irregularities, we have in Bennu a scenario where flattened surfaces coexist with micron-scale irregularities with different scattering behaviors, leading to another level of complexity for future photometric modeling for Bennu.

Phase Reddening
Phase reddening is a widely observed but poorly understood phenomenon for solar system objects, occurring on surfaces from basaltic, such as those of the Moon (e.g., Hapke et al. 2012), Mercury (e.g., Warell & Bergfors 2008, and Vesta (e.g., Li et al. 2013c;Schröder et al. 2014b); to ordinary chondritic, such as those of S-type asteroids (e.g., Clark et al. 2002;Kitazato et al. 2008); to carbonaceous, such as the surface of Ceres (Reddy et al. 2015;Ciarniello et al. 2017;Li et al. 2019;Longobardo et al. 2019) and Bennu; as well as the surfaces of comets (e.g., Fornasier et al. 2015;Feller et al. 2016) and icy moons (e.g., Buratti et al. 1990). It is also easily reproduced in the laboratory using either meteoritic samples or terrestrial analogs (e.g., Gradie & Veverka 1986;Beck et al. 2012;Potin et al. 2019). Despite its long recognition and some efforts in laboratory experimental studies (e.g., Beck et al. 2012;Schröder et al. 2014a;Pilorget et al. 2016), the exact physical origins of phase reddening are unclear. Some analyses have suggested that phase reddening could be caused by increasing contributions from multiple scattering with increasing phase angle for spectrally red surfaces (e.g., Hapke et al. 2012;Li et al. 2013c). But the existence of phase reddening on dark objects such as Ceres (Li et al. 2019) and Bennu, whose spectra are neutral and blue, respectively, and where the contribution of multiple scattering is negligible, is not consistent with this explanation. Laboratory experiments have suggested that phase reddening could originate from internal and external structures on the surfaces of scattering particles (Pilorget et al. 2016), or originate from microscale roughness that arises from fairy castle structures of loosely packed particles or surface roughness of particles (Schröder et al. 2014a). In any case, it appears that phase reddening is primarily controlled by the physical properties of regolith particles rather than composition (Beck et al. 2012), although the latter can contribute in cases of a diverse mineralogy (Jost et al. 2017).
We found that the phase-reddening map of Bennu does not seem to correlate with that of any other quantity (e.g., Figure 18(e)). Such a lack of correlation is consistent with phase reddening on Bennu's surface being controlled by microscale structure on the order of microns (Schröder et al. 2014a) but not by composition, space weathering, or macroscopic roughness (greater than tens of microns). If fine particles exist on Bennu, then the microscale roughness due to the formation of fairy castle structures (Schröder et al. 2014a) could be the cause of phase reddening. Otherwise, the phase reddening of Bennu's surface might be related to the microscopic roughness on the surface of particles at micron scales. If this is the case, then Bennu would be covered by particles that have rough surfaces at scales from microns to centimeters. Our findings are consistent with the results by Fornasier et al. (2020).

Latitudinal Trend and Space Weathering
On a global scale, an overall correlation exists among albedo, color, and thermal inertia, and a weak inverse correlation is possible between albedo and roughness (Figures 18(a)-(c)). The global correlation between high albedo and blue color has been discussed by DellaGiustina et al. (2020), who showed that both the primordial heterogeneity of boulders and space weathering of dark materials (regolith and boulders) contribute to such a correlation. Their work also indicated that particle size could affect the overall color trend of Bennu but may not dominate.
Another type of correlation is the latitudinal trends, as we will discuss in detail here. The equatorial region of Bennu appears to be relatively red in the spectral range of about 1-2 μm (Figures 11(d), 12(a), and 19), has higher thermal inertia (Rozitis et al. 2020b), and probably has slightly weaker phase reddening than mid-to-highlatitude areas (Figures 14(d), 15(a), and 19). The latitudinal color trend that we observed is consistent with previous analyses of OVIRS data (Barucci et al. 2020;Fornasier et al. 2020). In addition, Praet et al. (2021) reported a latitudinal trend in hydration, with low hydration in the equatorial region and generally increasing with latitude likely associated with preferential modification of Bennu's equator by space weathering, heating, and/or other processes.
As introduced in Section 1, ongoing latitudinal mass movement from higher latitudes toward the equator occurs on Bennu Jawin et al. 2020) as a result of continuous YORP spin-up ) and the geophysical dynamics on Bennu ). The equatorial region may also preferentially accumulate particles ejected from the global surface (McMahon et al. 2020). Therefore, the surface of Bennu continuously refreshes itself, with potentially more fresh material exposed at higher latitudes than in the equatorial region. The relatively red color of the equatorial belt in the NIR is consistent with a spectral reddening effect of space weathering, in the general sense of processes associated with longer or more intense exposure than higher-latitude regions, on Bennu in the 1-2 μm range.
However, if the redder equatorial region in the NIR is due to space weathering, why is such a red region lacking in the visible wavelengths, 0.44-0.89 μm (DellaGiustina et al. 2020; Figure 13(b))? Using principal component analysis of the MapCam colors, DellaGiustina et al. (2020) showed that albedo dominates the variations on Bennu's surface, while spectral slope variations are secondary. Their work suggested that compositional variations account for albedo differences on Bennu, consistent with findings by Sen et al. (2021), whereas space weathering-related spectral changes were a secondary factor. Thus, space weathering-related spectral effects in the visible wavelengths appear most prominently in local, compositionally uniform areas with different exposure ages.
Using OVIRS data, we confirmed that the equatorial region is relatively blue in the visible wavelengths (Figures 11(b), 12(a), and 13(a)), in line with results from MapCam (DellaGiustina et al. 2020). DellaGiustina et al. (2020) proposed that this enhanced bluing may arise from magnetite, which has an absorption feature at 0.55 μm, or graphitized carbon, which shows an upturn in the near-UV region. Indeed, Simon et al. (2020a) showed that the 0.55 μm band is enhanced along the equatorial region. The spectral matches to this band include magnetite, goethite, troilite, and certain phases of graphitized carbon, with magnetite being the best candidate (Simon et al. 2020a), which has been confirmed by OTES observations in MIR (Hamilton et al. 2019). Spectral modeling has also suggested the presence of magnetite and troilite on Bennu, although they do not seem to be concentrated along the equator (Trang et al. 2021). These species could all be related to space-weathering processes on carbonaceous chondrite materials and result in bluing in the visible wavelengths (e.g., Weisberg et al. 2006;Thompson et al. 2019).
If the latitudinal color trend on Bennu is due to space weathering, the colors of the equatorial region relative to the global average may reveal the differing spectral effects of weathering in distinct wavelength regions. In the NIR, weathering causes a strong reddening that dominates over composition-related color variations. In the visible, the spectral change is relatively weak, and the global color distribution is dominated by compositional heterogeneity. In the blue to visible wavelengths, spectral bluing is dominated by the enhancement of the 0.55 μm band, possibly due to the production of magnetite.
The explanation of relatively weak phase reddening in the equatorial region is uncertain, owing to the poor understanding of the phase-reddening effect and the uncertain abundance of fine particles on Bennu. Our best assessment is that it may imply some relationship between the micron-scale structure of Bennu's regolith particles and the latitudinal mass movement of Bennu's surface materials toward the equator.

North-South Asymmetry
Surface roughness and phase reddening both differ between the northern and southern hemispheres of Bennu. The southern hemisphere is rougher overall than the northern hemisphere (Figures 9(d) and (e)). Phase reddening is weaker in the southern hemisphere, with a magnitude similar to that of the equatorial band (Figures 14(a), 15(a), and 19). Despite the noise in the phase-reddening maps, the relationship appears to be stronger than the noise level. Fornasier et al. (2020) independently confirmed this north-south asymmetry with a different approach, where they directly measured the spectral slopes of the craters and rocks included in their study at various phase angles and then fitted the phase-reddening parameters.
Like the latitudinal trend in color and phase reddening that we discussed in the previous section, we attribute this northsouth asymmetry in roughness and phase reddening to the asymmetry in the surface geophysical environment and the resultant surface geology. Figure 18. Density plot of various parameter pairs and the corresponding correlation coefficients R showing their correlation or noncorrelation. Albedo is correlated (or anticorrelated) with thermal inertia (a), thermal roughness (b), and spectral slope (c) to some extent, whereas no correlation between albedo and phase slope (d) or phase reddening (e) is evident. Thermal roughness appears to be strongly correlated with phase slope (f).
The majority of large boulders (>30 m) on Bennu sit in the southern hemisphere Jawin et al. 2020). As a result, the latitudinal mass movement in the southern hemisphere is more obstructed by boulders, whereas the northern hemisphere shows more evidence of boulder dynamics and surface flow. This difference leads to an overall rounder shape for the southern hemisphere of Bennu, whereas the northern hemisphere is steeper with higher gravitational slope Daly et al. 2020). Considering that the equatorial band on Bennu also has a relatively low gravitational slope, there may be a correlation between gravitational slope and phase reddening, where higher gravitational slope, or potentially more mass movement, is correlated with stronger phase reddening.
If phase reddening is caused by fine particles forming fairy castle structures (Schröder et al. 2014a(Schröder et al. , 2014b, then perhaps a high gravitational slope favors the formation of such structures or a low gravitational slope favors the destruction of such structures by helping the particles settle. Schröder et al. (2014a) demonstrated that sprinkling, shaking, tapping, or pressing a sample all easily affected the magnitude of phase reddening that they measured. On the other hand, if phase reddening is controlled by the surface texture of particles or rocks at the micron scale, then perhaps more mass movement causes more collisions and friction between particles, creating a rougher texture. Or, weaker mass movement helps remove the microscale roughness by abrading the surface of the particles. Despite the highly uncertain physical interpretations of phase reddening, the correlation between phase reddening and gravitational slope is nonetheless plausible for Bennu's surface.
The north-south asymmetry in photometric roughness on Bennu might be correlated with the geological differences in terms of boulders and mass movement. If, as we discussed before, the roughness on Bennu is similar at scales greater than tens of microns, then the concentrated distribution of large boulders in the southern hemisphere, and the concentration of small rocks that might be broken down from those large boulders, could be an explanation for the higher roughness in the southern hemisphere.

Comparison with Other C-complex Asteroids
The bulk photometric properties of Bennu are generally similar to those of other primitive asteroids, in particular (253) Mathilde (Clark et al. 1999) and (162173) Ryugu (Tatsumi et al. 2020), which have similar albedo and phase function. Cometary nuclei (Li et al. 2007a(Li et al. , 2007b(Li et al. , 2009(Li et al. , 2013a(Li et al. , 2013bFornasier et al. 2015) also have an albedo similar to Bennu of about 0.04, but their phase functions are slightly steeper, presumably related to their more porous surfaces causing more prominent shadowing than on Bennu.
Most relevant to Bennu is Ryugu, the target of JAXA's Hayabusa2 mission (Watanabe et al. 2017), which is also covered by aqueously altered carbonaceous materials, as indicated by the widespread 2.7 μm spectral absorption on its surface (Kitazato et al. 2019). Its top-like shape is also similar to Bennu's, indicating the operation of the YORP effect (Watanabe et al. 2019).
The overall photometric properties of Bennu are similar to those of Ryugu, with a slightly higher albedo in the visible wavelengths and a slightly shallower phase function Figure 19. A 3D view of the spectrophotometric maps. The maps are overlaid on the base map of Bennu (Bennett et al. 2020) and visualized in the Small Body Mapping Tool (Ernst et al. 2018). The left column is the view from the -x direction with a central longitude of 0°, and the right column is from the +x direction with a central longitude of 180°. The normal albedo and phase slope maps are at 0.55 μm. The spectral slope and phase-reddening maps are for the spectral range 0.5-2.0 μm. The ranges of the color bars are the same as those in Figures 9, 11, and 13. Normal albedo is a dimensionless quantity, and the units of other maps are mag deg −1 for phase slope, μm −1 for spectral slope, and μm −1 deg −1 for phase reddening.  (Tatsumi et al. 2020). Phase reddening of Bennu in the visible wavelength (Figures 13 and 14) is weaker than that of Ryugu (Tatsumi et al. 2020), as also reported by Fornasier et al. (2020) and Zou et al. (2021). Equatorial banding also exists on Ryugu, although in a different way than is observed on Bennu. The equatorial region of Ryugu appears to be bright and blue in the visible wavelengths (Sugita et al. 2019;Tatsumi et al. 2020) but shows no significant variation from other areas in the NIR (Kitazato et al. 2019;Pilorget et al. 2021) and has relatively low thermal inertia and roughness (Morota et al. 2020;Shimaki et al. 2020). Evidently, such equatorial banding is also associated with mass movement on Ryugu's surface (Sugita et al. 2019).
The geophysical conditions and resultant mass movement on Ryugu are different from those on Bennu in two respects. First, the global average gravitational slope on the surface of Ryugu is about 12° (Watanabe et al. 2019), the same as the average slope on Bennu within the equatorial band but much lower than that found in other areas on Bennu (18°; Scheeres et al. 2019). Therefore, the correlation on Bennu's surface between a relatively low gravitational slope and a low level of phase reddening (Sections 4.4) does not apply to the comparison between Bennu and Ryugu. This might indicate that the correlation we find on Bennu might be associated with specific physical or mineralogical properties and cannot be generalized to other asteroids.
Second, the equatorial region of Ryugu represents a gravitational high. Correspondingly, the direction of mass movement on Ryugu is from the equator toward high latitudes (Sugita et al. 2019), the opposite direction of mass movement on Bennu. Therefore, freshly exposed underlying material dominates the equatorial region of Ryugu, and it is brighter and bluer than the older surface materials that have been reddened by solar heating or space weathering (Morota et al. 2020). This interpretation is also consistent with the relatively blue color of small, young craters on Ryugu. Therefore, we might be able to generalize the equatorial banding in the spectral color as evident on both Bennu and Ryugu to other similarly sized and shaped near-Earth asteroids, although the details may differ.

Summary
Our global photometric modeling confirms that Bennu exhibits the typical overall photometric characteristics of a dark, primitive solar system object. At 0.55 μm, the geometric albedo is 0.046, and the Bond albedo is 0.020, according to the Akimov empirical model with an exponential phase function, using the OVIRS data at phase angles 7°-100°. These values are higher than previous ground-based results by about 10% (Hergenrother et al. 2013;Takir et al. 2015), which is attributable to the absolute flux calibration of OVIRS. The uncertainty in the albedo that we derive is 15%. The slope of Bennu's phase function is 0.024 mag deg −1 for the surface phase function, translating to 0.038 mag deg −1 for the diskintegrated phase function, with an uncertainty of 30%. The phase slope of Bennu that we measure is slightly lower (shallower) than in previous ground-based observations (Hergenrother et al. 2013;Takir et al. 2015) and the OCAMS Approach-phase data ) but consistent within the uncertainties. Bennu has a global average phasereddening parameter of 4.3 × 10 −4 (μm deg) −1 , and its phase reddening decreases with wavelength. Bennu's phase reddening is weaker than that of Ryugu and other carbonaceous asteroids, such as Ceres. Bennu's surface is spectrally blue with a global average spectral slope of −0.0030 μm −1 in 0.5-2.0 μm or −0.7% per 100 nm when normalized to the middle of this spectral range.
Bennu's surface is covered by a well-mixed regolith with rocks of heterogeneous compositions and properties and a wide range of albedos at all length scales, which has not been observed on other asteroids, including Ryugu. Such a regolith, combined with the complicated local topography caused by rocks, made it difficult to reliably constrain the Hapke roughness parameter. The high uncertainty in the modeled Hapke roughness causes high uncertainty in the modeled single-particle phase function, which in turn affects the reliable determination of the SSA. The empirical model with Akimov disk function and an exponential surface phase function seems to fit Bennu's global photometric data, although the residual is higher than in previous modeling of other asteroids, including Ceres and Vesta.
Regional spectrophotometric mapping with the Akimov empirical model yields normal albedo and phase slope maps of Bennu at all spectral channels of the OVIRS instrument, although the fitting quality at wavelengths >2.7 μm is potentially affected by thermal tail removal. The albedo map that we derived at 0.55 μm is consistent with those previously derived from OCAMS (Golish et al. 2021b) and OVIRS (Zou et al. 2021) data. Dark areas are generally correlated with large dark boulders and boulder-rich areas in the southern hemisphere. No global correlation appears between albedo and phase slope, although some correlation is evident at regional scales, suggesting that roughness variations likely dominate the phase slope variations.
Two types of correlations between albedo and spectral slope exist on Bennu. On a global scale, high albedo generally correlates with a relatively blue spectral slope, as well as high thermal inertia. This correlation could be attributed to the primordial heterogeneity of rocks on Bennu and the associated albedo, color, and thermal inertia differences. On the other hand, our analysis revealed a latitudinal trend in spectral slope and phase reddening on Bennu. The equatorial region is associated with the reddest (least blue) spectral slope in the 0.5-2.0 μm range, as well as relatively weak phase reddening. The spectral slope decreases toward higher latitudes. Phase reddening also shows a slight north-south asymmetry, where the northern hemisphere has a slightly stronger reddening than the south.
We interpret the latitudinal trend of spectral slope and phase reddening with respect to the geophysical conditions and consequent geology on Bennu. The continuous YORP spin-up of Bennu produced gravitational highs in the polar and midlatitude regions and gravitational lows in the equatorial region . The surface material on Bennu flows toward the equator (Jawin et al. 2020). Ejected particles are also more likely to collect in the equatorial region (McMahon et al. 2020) owing the structure of Bennu's rotational Roche lobe . Therefore, the equatorial region collects geologically old materials, and higher-latitude regions contain progressively more abundant young, freshly exposed materials. Therefore, the latitudinal trend of spectral slope in the 1-2 μm region that we observed is possibly a result of space weathering. Comparison with the color maps in the 0.55-0.89 and 0.44-0.55 μm region (DellaGiustina et al. 2020) suggests that the spectral effects of space weathering on Bennu vary with wavelength. The equatorial band of weak phase reddening and its slight northsouth asymmetry could also be associated with the direction and predominance of surface mass flow and the north-south asymmetry of the geophysical conditions on Bennu ), but the exact cause is unclear, owing to the poor general understanding of phase reddening on planetary surfaces, especially for those with dark surfaces.