Regional Photometric Modeling of Asteroid (101955) Bennu

We present a regional photometric analysis of asteroid (101955) Bennu, using image data from the MapCam color imager of the OSIRIS-REx Camera Suite (OCAMS). This analysis follows the previously reported global photometric analysis of Bennu, which found that Bennu’s roughness was difficult to photometrically model owing to unresolved surface variation. Here we find that, even with a high-resolution shape model (20 cm per facet) and automatic image registration (<1 pixel error), Bennu remains a challenging surface to photometrically model: neither a suite of empirical photometric models nor the physically motivated Hapke model were able to eliminate the scatter in the data due to pixel-scale variations. Nonetheless, the models improved on the global analysis by identifying regional variations in Bennu’s photometric response. A linear empirical model, when compared with independent measures of surface roughness and albedo, revealed correlations between those characteristics and phase slope. A regional Hapke analysis showed the same structure in its single-scattering albedo and asymmetry factors; although the Hapke parameters were loosely constrained, complicating interpretation of their spatial variation, the regional variation in relative parameter sensitivity also correlated with shallower phase slope, higher albedo, and less macroscopic roughness.


Introduction
The Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer (OSIRIS-REx) mission studied the surface of Bennu for 2 yr before collecting a sample from a site called Nightingale in Hokioi Crater on 2020 October 20. OSIRIS-REx spacecraft data confirmed ground-based observations of Bennu as a low-albedo, B-type asteroid with a slightly blue spectrum (Clark et al. 2011;Hergenrother et al. 2013;Lauretta et al. 2019;DellaGiustina et al. 2020). A global photometric analysis (Golish et al. 2021a) fit Bennu's photometric response to a suite of empirical models using color images acquired by MapCam of the OSIRIS-REx Camera Suite (OCAMS; Rizk et al. 2018). That work also verified ground-based expectations (Hergenrother et al. 2013;Takir et al. 2015) of a relatively steep phase slope commensurate with Bennu's low albedo and highly shadowed surface. It also found that a Robotic Lunar Observatory (ROLO; Buratti et al. 2011) formulation was required to capture the subtle opposition effect, and that measurable phase reddening occurs across the ∼400 nm wavelength range spanned by MapCam's four color bands (Golish et al. 2021a). Golish et al. (2021a) treated Bennu as a homogenous surface, using the median reflectance and photometric angles of every image (each of which covered ∼3/4 of Bennu's hemisphere) when fitting its photometric models. However, Bennu has a heterogeneous surface with notable variations in surface texture (Walsh et al. 2019;Daly et al. 2020;Bennett et al. 2021), color (DellaGiustina et al. 2020, albedo (DellaGiustina et al. Golish et al. 2021c), and composition (DellaGiustina et al. 2020; Rozitis et al. 2020;Simon et al. 2020). Moreover, the surface has experienced ongoing changes due to particle ejections, mass movement, impacts, and thermal fracturing (e.g., Ballouz et al. 2020;Hergenrother et al. 2020 and references therein;Jawin et al. 2020;Molaro et al. 2020).
In addition, Golish et al. (2021a) demonstrated that Bennu's rough surface introduces significant scatter to the reflectance and photometric data used to model it. The photometric response of a surface is determined by the phase angle (α), incidence angle (i), and emission angle (e) of the terrain in a given observation geometry. Some models also use photometric latitude (l) and longitude (b), which are described in Kreslavsky et al. (2000): As such, uncertainty in the photometric angles used by the model directly affects the accuracy of the model. Golish et al. (2021a) suggested that higher-resolution shape models, which provide a more accurate measure of the photometric angles, would be required to understand the photometric response of the surface.
In this work, we address Bennu's heterogeneity and roughness by analyzing its surface regionally using a shape model with increased resolution (20 cm versus. the previous 80 cm per facet; Daly et al. 2020). We fit the suite of empirical models previously used by Golish et al. (2021a) and perform a regional Hapke analysis of the surface. The spatial variation of Bennu's photometric properties (phase slope, phase reddening, and albedo) and their correlation with other characteristics of Bennu's surface provide a more detailed understanding of the photometric behavior of this heterogeneous body.

Observations
Image data used in this study were acquired by MapCam, the medium-angle color OCAMS imager (Rizk et al. 2018). MapCam imaged Bennu's surface in four narrowband (0.060-0.090 μm wide) color filters with center wavelengths at 0.473, 0.550, 0.698, and 0.847 μm. We also used panchromatic image data acquired with MapCam's broadband filter (∼0.3 μm wide, centered at 0.646 μm). The OCAMS calibration pipeline (Golish et al. 2020) corrected these data for known sources of error and converted them to units of I/F, also referred to as radiance factor (RADF; Li et al. 2015). I is the scattered radiance (W m −2 μm −1 sr −1 ), and F is the radiance from a Lambertian surface normally illuminated by the Sun.
We used the same image set as the global photometric analysis (Golish et al. 2021a) with a few modifications. The majority of the photometric data were from the Equatorial Stations observation campaign (DellaGiustina et al. 2018;Lauretta et al. 2021), which acquired images at phase angles ranging from 7°to  130 , with a pixel scale of ∼33 cm pixel -1 ( Table 1). The spacecraft was positioned over Bennu's equator for every Equatorial Station. We refer to each Equatorial Station by its local solar time (LST); the 3 pm station, for example, corresponds to phase angles of ∼45°. Golish et al. (2021a) did not include the high-phase-angle (∼130°) data owing to excess noise. In this analysis, with a more accurate shape model and registration to that model, we did include the high-phase-angle data. As in the global analysis, we supplemented these data with images from the Preliminary Survey mission campaign. These images were acquired from farther away (pixel scale of ∼1 m pixel -1 ) and did not provide global coverage, but they did observe phase angles (38°-90°) that were not imaged during Equatorial Stations (Table 1).
Unlike the global photometric analysis, we also included color MapCam images from flybys (FBs) 2(b), 7(a), and 7(b) of the Baseball Diamond campaign (DellaGiustina et al. 2018;Lauretta et al. 2021). MapCam acquired the FB2b images at ∼7°phase with an equatorial subspacecraft position, similar to those from the 12:30 pm Equatorial Station. MapCam acquired the FB7a and FB7b images, however, with midlatitude subspacecraft positions. These images were the only global color observations acquired by OSIRIS-REx that imaged the midlatitudes at low emission angles, providing unique photometric views of the surface.
The low-phase-angle data collected during the Approach mission campaign (Lauretta et al. 2021), which were included in the previous global photometric analysis, were not included here because Bennu subtended a small number of pixels (<50) in those images. Therefore, the lowest phase angles included in this analysis are ∼7°. Accordingly, we did not attempt to model Bennu's opposition surge regionally. Instead, we took advantage of previous analyses and used the opposition surge terms derived therein for both empirical (Golish et al. 2021a) and Hapke (Hergenrother et al. 2019) modeling.
We used >4000 MapCam images from the panchromatic and four color bands (Table 1). Figure 1 depicts maps of the minimum and maximum photometric angles, measured in 1°×1°latitude/ longitude bins with the v band; all five filters had roughly equivalent coverage. The data set has excellent phase angle coverage over the majority of the surface-only phase angles under 7°are unrepresented. High phase angles are sporadically missing across the surface owing to Bennu's rugged terrain; at high phase angles, the Sun does not illuminate large swaths of the surface. Thus, some regions were observed only at 90°phase. The polar regions were typically limited to 45°phase. Incidence angles were determined by Bennu's terrain and orientation. Bennu has a very low obliquity (∼177°.6; Lauretta et al. 2019), so the incidence angles were lowest at the equator. Minimum incidence angles increased sharply in the regions past the edge of the equatorial bulge . MapCam rarely imaged the surface with emission angles below 10°, again owing to Bennu's terrain. However, it did image the majority of the surface with emission angles down to  20 . The maximum observed emission angles, up to ∼65°, were likewise determined by the terrain. Histograms of these data ( Figure 2) further illustrate the fraction of the surface imaged with a given minimum or maximum photometric angle.
One factor not illustrated by these plots is the diversity of emission and incidence angles within a phase angle range. As noted above, we only acquired low-emission views of the midlatitudes in FBs 7(a) and 7(b), which had phase angles between 31°and 45°. So although the maps show phase coverage from 7°to 135°and emission angle coverage from 20°to 70°, there were many portions of the surface where those are not both true at the same time. Emission (and incidence) angle coverage is primarily important for understanding the surface's disk function. The disk function has a weak or nonexistent phase relationship, depending on the formulation. Conversely, phase angle coverage is important for modeling the surface's phase response. As such, we expect that this limitation does not dramatically affect our results.

Automatic Image Registration
The previous global photometric analysis demonstrated that shape model resolution has a direct effect on the quality of photometric modeling and correction (Golish et al. 2021a). This behavior is not unique to Bennu; to model the dependence on phase, incidence, and emission angles, we must know those angles as accurately as possible. What is particular about Bennu, and presumably other rubble-pile asteroids such as Ryugu (Tatsumi et al. 2020), is the extraordinary roughness that produces high spatial frequency variations in emission and incidence angles. Registering the images to a high-resolution shape model is required so that the model accurately represents the variation seen in the images.
Fortunately, Bennu's surface has been measured with the OSIRIS-REx Laser Altimeter (OLA; Daly et al. 2017), which made it possible to construct the highest-resolution shape model of any planetary body (Seabrook et al. 2019;Barnouin et al. 2020;Daly et al. 2020). This model, referred to as OLA v20, has a ground sample distance (GSD) of 5 cm at its native resolution. The majority of the image data used in this study had pixel scales ranging between 25 and 33 cm pixel -1 . Computational limitations did not allow us to take advantage of the full 5 cm shape model, which has over ∼1.5 billion facets. However, the authors of Daly et al. (2020) also resampled the OLA v20 shape model into tiles with 20 cm GSD; we adopt this version here.
Each of the 20 cm GSD tiles covers ∼15°×15°of latitude/longitude and is composed of >600,000 facets; the complete global set is~150 million facets. The ray-tracing time associated with this number of facets for >4000 images was computationally prohibitive. Therefore, we used a custom Python script to automatically identify all shape model tiles represented in a given image (e.g., not out of field or on the backside of the asteroid). Approximately 50 20 cm tiles were required to cover a typical Equatorial Stations image (example in Figure 3). Though handling this many facets (>30 million) was complex, it was feasible with cluster computing resources.
We registered the 20 cm tiles to the images in the U.S. Geological Survey's Integrated Software for Imagers and Spectrometers 3 (ISIS3; Keszthelyi et al. 2013), with which we calculated the photometric angles of each pixel in each image. The initial alignment of the images to the shape model was provided in the a priori SPICE kernels (Acton et al. 2018) developed by the OSIRIS-REx navigation team. However, these kernels were only accurate to 1-2 m in ground distance, which corresponds to an offset of 3-6 pixels in an Equatorial Stations image. As such, our first processing step was to align the images to their shape model tiles.
We measured the initial offset between an image and the shape model by creating simulated images (DellaGiustina et al. 2018). In ISIS3, using the a priori kernels, we calculated the phase, incidence, and emission backplanes for a given image. From those backplanes, we created a photometric simulation of the image using the Lommel-Seeliger photometric model presented in Golish et al. (2021a). The primary goal of the photometric simulation was to mimic the disk function variations across the image, which can be dramatic given Bennu's rough surface, though the simulations also accommodated gentle variations in phase angle across the image. Importantly, we also used ISIS3ʼs ray-tracing engine to determine the regions of the surface that did not have a direct  line of sight to the Sun (shadows), which we masked out of the simulated image. At all phase angles, shadows were critical for registering the images. The accuracy of the shape model, particularly the shape and height of surface features , directly determined the fidelity of the shadows. A representative image, acquired at 2019 April 25 17:36:46 UT during the 3 pm Equatorial Station, is shown in Figure 4(a); its simulated counterpart, created using a priori kernels, is shown in Figure 4(b). An alignment mark is overlaid on both images to indicate the ∼6 pixel offset between the two images.
Typically, we would update the registration of a set of images to a shape model using the photogrammetric techniques described in Bennett et al. (2021) and Edmundson et al. (2020). Those processes used automatic feature identification to generate thousands of tie points to produce a control network that tied all images to each other. However, these automated schemes were most effective for images at similar observation conditions. Automatically identifying a given feature in images taken at different local times of day is much more challenging, especially because of Bennu's rough terrain. Hence, we developed a new automatic method to register the images used in this study.
Our automatic registration procedure used the same photometric simulations described above as a photogrammetric reference with a series of ISIS3 applications. This process began with the findfeatures ISIS3 application, which used a Features from Accelerated Segment Test (FAST; Rosten & Drummond 2005) detector and a Scale Invariant Feature Transform (SIFT; Lowe 2004) extractor to identify a control network of tie points between the actual and simulated images. The locations of these tie points were then refined using pointreg, an intensity-based pattern-matching ISIS3 application. We configured pointreg uniquely for each phase angle grouping, as the quantity of shadows strongly affects the number of valid pixels in the search window. Finally, we used the bundle adjustment routine jigsaw to solve for the updated pointing that minimized the error in the bundle of rays described by the control network (Edmundson et al. 2020). After updating the pointing and recalculating a simulated image, the image matched its simulation to better than 1 pixel (Figure 4(c)).
A disadvantage to this method is that we solved for the registration of each image independently. There was no imageto-image network such as is typically preferred for mosaicking, when images are projected to a shared coordinate system. However, for the purposes of generating photometric backplanes, image-to-image alignment was not necessary. Only registration to the shape model affects photometric modeling. An advantage of this method was that we could automatically apply it to all images in our data set. The automatic registration was equally effective for low-phase-angle images with minimal shadows (Figures 4(d)-(f)) and high-phase-angle images with extreme shadows (Figures 4(g)-(i)). This automatic registration method is broadly applicable to images of any planetary surface where a sufficiently high-resolution shape model (∼pixel scale) is available.

Latitude and Longitude Binning
After registering all images to the shape model, we recalculated the photometric backplanes for each image. We calculated geographic latitude and longitude backplanes, based on the location of the rays' intersection with the shape model. We used these backplanes to divide the images into latitude/ longitude bins, calculating the median reflectance, phase, incidence, and emission for each bin. We calculated photometric latitude and longitude for each bin based on Equations (1) and (2). Because the photometric models are more sensitive at high emission and incidence angles and the shape model is less accurate on steep slopes , we ignored pixels observed with emission or incidence angles greater than  82 . To mask out shadows and increase the signal-to-noise ratio of the data, we also ignored pixels with I/F values less than 0.0001.
We merged bins from all images and required that each latitude/longitude bin contain at least 15 data points (i.e., data from 15 different images) and that each bin had complete phase angle coverage between 7°and 90°. This removed portions of the surface with insufficient diversity in phase angle ( Figure 1). Based on the results of this data binning, we adopted the 1°×1°latitude/longitude grid depicted in Figure 5. Selecting a smaller bin size reduced the photometric coverage. In particular, the number of bins that maintained a maximum phase over 90°was reduced and the number of bins with low minimum incidence angles and high maximum emission angles decreased. We performed regional photometric modeling on individual bins.

Photometric Models
The previous global photometric analysis of Bennu fit empirical models to the data (Golish et al. 2021a). These models were composed of pairs of phase functions (e.g., Exponential and ROLO) and disk functions (e.g., Lommel-Seeliger and Lunar-Lambert). In Golish et al. (2021a), the models were referred to by both their disk function (e.g., Lommel-Seeliger) and their phase function (e.g., ROLO). Though the originators of each model chose the phase and disk functions in tandem, we have the flexibility to investigate them independently. Therefore, in this work we broadened that analysis to include all combinations of the phase and disk functions investigated in Golish et al. (2021a), which we list in Tables 2 and 3.

Global Disk Function
Three of the four disk functions considered in this study are parameterized; Lommel-Seeliger is fixed. Golish et al. (2021a) found significant scatter in the disk function solutions, even at low phase. The OLA v20 20 cm GSD shape model tiles used in this work represent a4 increase in resolution and a more accurate measure of surface roughness. As such, before performing a regional analysis, we repeated the global disk function analysis to evaluate the improvement from a higherresolution shape model.
Variation in the normal albedo across Bennu's surface is substantial (Golish et al. 2021c) and can account for some scatter in the data. To minimize this effect, we created albedo simulants of every MapCam image. We used the backplanes described in Section 2.2 and the normal albedo map from Golish et al. (2021c) to assign a normal albedo to each pixel in an image (both data sets were registered to the OLA v20 shape model). We then normalized the albedo image and divided it into the corresponding MapCam image. Finally, we removed phase angle variation across the image by applying the corresponding phase function correction from the previous global analysis (Golish et al. 2021a). After this series of corrections, each image theoretically represented disk function variation only.
For direct comparison, we repeated the techniques performed in Golish et al. (2021a) for these albedo-and phase-corrected images. That is, we reformulated the Lunar-Lambert (Equation (3)), Minnaert (Equation (4)), and Akimov (Equations (5)- (6)) disk functions such that their free parameters (L(α), k(α), and η(α), respectively) were the slope of a linear relationship dependent on the photometric angles (i,e,b,l) and albedo-and phase-corrected reflectance d(i,e)). These formulations are not fundamentally different from solving the disk functions directly, but they allowed for easier visualization: Despite the higher-resolution shape model, more accurate image registration, and albedo/phase removal, we still find significant scatter in the disk modeling data. Figure increasing scatter with increasing phase angle. Only the Minnaert function maintains a positive R 2 value at all phase angles, implying some linear correlation. The parameters of these fits are uncertain. Therefore, we combined the solutions from every image in our data set to reduce the noise and to determine a phase relationship ( Figure 8). Unfortunately, the majority of images fit the Lunar-Lambert parameter L to 1, such that there was effectively no phase dependence. This was made clear by the examples in Figure 7: as the phase angle increased, the quantity of data acquired at low incidence and emission angles (the lower left corner of the plot) vanished. Correspondingly, the density of points along the left edge of the scatterplot increased and thus influenced the solution more strongly than the rest of the data. The Akimov disk function, for the reasons discussed above, did not produce a meaningful relationship. The Minnaert disk function was the only one with a positive R 2 value and a phase dependence that is plausible.
Each of the three different empirical disk functions that we attempted to fit proved insufficient, due to the very noisy data set. To better visualize this, we plotted the incidence and emission angles from the 12:30 pm LST image against the measured I/F ( Figure 9) and overlaid a three-dimensional surface representing the Lommel-Seeliger disk function as a reference. At every combination of incidence and emission, there are multiple measurements over a range of I/F values. The standard deviation of the I/F variation around the Lommel-Seeliger reference (0.0038; ∼10% of the median) represents the noise in the data. We visualize this noise by calculating the histogram of the image with and without removing a Lommel-Seeliger disk function (Figure 10(a)). Zooming in on a 50×50 pixel portion of the surface ( Figure 10(b)) shows that even after Lommel-Seeliger correction was applied the majority of the high-frequency variation remained. The ability of the model to describe the surface worsens as phase angle increases, where more shadows increase the I/F variation. Moreover, images acquired at higher phase angle necessarily limit themselves to higher incidence and/or emission angles, which increase the scatter in the data. For example, at 3 pm LST (Figure 11), the standard deviation increases to 0.0085 (∼20% of the median). These results illustrate that the higher-resolution shape model does not eliminate the noise in our photometric data.

Regional Disk Function
A possible cause of our inability to fit a disk function to these images might be regional variation in that disk function. In the previous section, we attempted to fit a single function to Bennu's surface. As discussed, Bennu has a heterogeneous surface, and it is reasonable to consider a spatially variant disk function. We used the latitude/longitude bins described in Section 2.3 to analyze the disk function regionally. That is, we used every image that included a given latitude/longitude bin, at all measured phase angles. For each bin, we solved for the free parameters of the Lunar-Lambert, Minnaert, and Akimov disk functions, including the phase angle dependence. The data were still sufficiently noisy that allowing the Lunar-Lambert partition function L(α) to have a third-order polynomial, as Figure 6. A representative 12:30 pm LST (7.5°phase angle) image (panel (a)), when corrected for albedo (panel (b)) and phase angle variation, should represent only disk function variation (panel (c)). However, removing albedo and phase narrows the I/F variation only slightly (panel (d)). Disk function parameter fits for the Lunar-Lambert (e) and Minnaert (f) models, where warmer colors indicate a higher density of data, are noisy but feasible. The fit for the Akimov model (g) is essentially unconstrained.
formulated by McEwen (1986), tended to allow L functions that were nonmonotonic. As such, we reduced the polynomial to first order (i.e., forced the Lunar-Lambert ξ and η to be 0).
We corrected for phase angle variation across the image using the phase function parameters published in Golish et al. (2021a). We did not correct for albedo, as there is much less  albedo variation across 1°×1°bins. We used a least-squares solver (MATLAB's lsqcurvefit) to optimize the free parameters for the (i,j)th bin (C i,j ) such that it minimized the difference between the modeled and measured reflectance (I F i j , ). The albedo of every bin (A i,j ) was also an output of this process, as an independent scalar multiplied by disk function:   . 3D plot of the incidence and emission angles from the 3 pm LST image against the measured I/F demonstrates higher noise at a larger phase angle. It also shows that low emission and incidence angles are not available at this phase angle, but a broader set of high photometric angles are.
We calculated the standard deviation of the difference between the disk function model and the image data for every latitude/ longitude bin using both the regionally solved parameters and the global parameters from Golish et al. (2021a), allowing for a direct assessment of the efficacy of regional modeling. Figure 12 depicts the disk function maps that resulted from this analysis. The top row depicts maps of albedo that, for the Lunar-Lambert and Akimov models, qualitatively match the published albedo map of Bennu (Figure 12(i); Golish et al. 2021c) The Lunar-Lambert albedo map has less contrast than the Akimov map, but it also has less high-frequency noise (individual bins that solved poorly). The Minnaert albedo map deviates noticeably from Bennu's known albedo and seems to be sensitive to topography.
The parameter maps (Figure 12, middle and bottom rows) depict interesting structure, particularly for the Akimov disk function. The Akimov parameter η tends to be higher in regions with rough, low-albedo boulders. The Akimov parameter q( = ha p a -q ), i.e., the exponent of the cosine of the photometric longitude (Table 3), is referred to as the smoothness parameter (Korokhin et al. 2007;Shkuratov et al. 2011) and has been shown to increase with increasing surface roughness on the Moon (Kreslavsky et al. 2000;Velikodsky et al. 2000Velikodsky et al. , 2002. Moreover, the smoother regions of the Moon (e.g., the maria) are also darker than the rough highland regions. That η is higher in rough (but dark) regions on Bennu is our first line of evidence that roughness is driving the photometric behavior of Bennu's surface. On the other hand, these data were corrected for phase angle using the global photometric solution from Golish et al. (2021a). These solutions did not include any regional phase function variation and therefore would not accommodate the steeper phase slope in darker and rougher regions of the asteroid. Therefore, what appears as η variation may simply be uncorrected phase variation. This does not invalidate the fact that regions of high η correlate with rough and dark features. However, it does expose an inherent degeneracy in the interpretation.
The Lunar-Lambert partition function, L(α), is nearly 100% Lommel-Seeliger (i.e., a = L 1 ( ) ) at low phase but tends to be lower (more Lambertian) at high phase in those same dark, rough regions of Bennu. This appears contrary to our expectation that darker surfaces are naturally single scattering and therefore tend to follow the Lommel-Seeliger law (Fairbairn 2005). However, the majority of Bennu's surface is among the darkest in the solar system Golish et al. 2021c); any Lunar-Lambert spatial structure is unlikely to be driven by albedo. Macroscopic roughness may itself induce multiple scattering, but we would still expect it to be minimal given Bennu's low albedo. However, as with the Akimov model's η term, uncorrected phase function variation is likely to manifest as variation in L(α) (Figures 12(e) and (h)).
Finally, the Minnaert b and k 0 maps have some structure that correlates with Bennu's albedo/roughness distribution (e.g., elevated b values in dark rough areas). However, their most prominent feature is a visible equatorial band. This equatorial band is likely an artifact, given that the Minnaert albedo map, which correlates poorly with Bennu's normal albedo map (Section 3.4), has the same band. Figure 13 plots maps of the standard deviation of the modeled versus measured reflectance using global (top row) and regional (bottom row) disk functions. The regional disk functions reduced the standard deviation (i.e., the model matches the data better) for all three models. We display the median standard deviation on the top of each map, which  (c)) qualitatively match Bennu's normal albedo map (panel (i); Golish et al. 2021c). The Minneart parameter maps (panels (d) and (g)) are dominated by an equatorial artifact, but the Lunar-Lambert (panels (e) and (h)) and Akimov (panel (f)) maps have structure that correlates with roughness and albedo trends on Bennu.
shows a global improvement of ∼20%-50%. Here we show the results for MapCam's v band. The other filters have similar behavior, though the median standard deviation varies by 15%-20%. Figure 14 reveals the mechanisms behind the improvement from global to regional; it plots the difference between the measured I/F (corrected for phase angle variation) and the I/F predicted by the global and regional disk function models for a representative latitude/longitude bin. The regional model data are clustered around zero because each bin has its own albedo value, whereas the global model has only a single value for the global albedo. The majority of the improvement in regional modeling was from this factor. However, the spread within each phase angle cluster (corresponding to data from different stations) is smaller, suggesting a better fit with the regional solutions. Given the improved fit and the spatially correlated parameter maps in Figure 12, the regional models are distinct improvements over their global counterparts.

Phase Function
We performed regional phase function modeling on the same 1°×1°binned data set used for regional disk function modeling. We corrected the measured I/F for disk function variation first and therefore modeled the equigonal albedo. We also corrected the data with all available disk functions (Table 3), including regional versions, and fit them against our four phase functions ( Table 2). For example, rather than forcing the Minnaert phase function to model the data corrected by the Minnaert disk function, the Minnaert phase function was fit to the data corrected by every disk function. For each latitude/longitude bin, we again fit the data using a leastsquares solver. Because we did not include Approach data in the regional analysis (owing to lack of resolved surface coverage; Section 2.1), we could not fit the ROLO opposition surge. Instead, we used the filter-specific opposition surge terms (C 0 and C 1 ) from the previous global analysis (Golish et al. 2021a). We solved phase functions for MapCam's five filters independently. The results are generally similar to each other; we primarily show the v filter (0.55 μm) results in the following sections.
Repeating this process over the entire surface produced a map of phase function parameters for each photometric model, including the albedo. Figure 15 plots the albedo term(s) for each phase function when solved with the Lommel-Seeliger disk function. Each map is paired with a plot correlating the albedos modeled by this analysis with the albedos measured in the normal albedo map (Figure 12(i); Golish et al. 2021c). The normal albedo map was created using PolyCam rather than MapCam data; we used it here as a reference to measure the efficacy of photometric solutions. Though the data should theoretically follow a 1:1 line, because the radiometric calibration of the two cameras is slightly different (Golish et al. 2020), we allowed the fit to have a slope not equal to 1. Figure 13. Maps of the standard deviation of the measured vs. modeled 0.55 μm data have lower values (i.e., better fits) when using regional solutions for all three disk functions. Figure 14. The difference between the measured and modeled I/F, at 0.55 μm, from the global (red asterisks) and regional (blue asterisks) disk functions for the bin at (29.5°, 179.5°) demonstrates the improvement offered by the regional solution.
The Minnaert phase function shows evidence of terrain, suggesting that observational biases are driving the solution.
Correspondingly, its correlation plot shows the poorest fit to the normal albedo map. The exponential phase function tends to return brighter albedos than the other models. The ROLO model's average albedo is better, but also with more scatter, leading to less correlation with the normal albedo map. This may be a result of the lack of a regional opposition surge term when modeling with ROLO. Each bin used Bennu's average opposition surge parameters. This is a limitation of ROLO without low-phase, high-resolution data. The albedo predicted by the linear model has the most internal consistency and the best correlation with the normal albedo map (Golish et al. 2021c).
Similarly, Figure 16 plots albedo maps for the linear phase function when solved with each disk function. Here the regional Minnaert disk function overestimates the albedo, particularly off the equator, perhaps driven by nadir imaging conditions at equatorial latitudes. The global Minnaert solution (not shown) overemphasizes the equator. This reinforces our conclusion from Section 3.3 that the Minnaert solution is not valid. The remaining disk functions have a similar qualitative appearance and quantitative correlation. The global versions of the Lunar-Lambert and Akimov functions (not shown) are similar to the regional versions.
The least-squares solver also calculated the squared norm of the residual between the measured and modeled data. We calculated the median of the squared residual for each phase/ disk function combination and compiled the results in Table 4. In line with our conclusions from the regional albedo models, the exponential phase function consistently performed ∼2× worse than other phase functions, regardless of disk function. Though the Minnaert phase function had residuals similar to the other phase functions, its poor correlation with the normal albedo of Bennu led us to discard that solution. The residuals are very similar for the ROLO and Akimov phase functions, with ROLO generally being slightly (but statistically insignificantly) better.
The global and regional Minnaert disk functions did not fit the data well, in line with the lack of correspondence between the Minnaert albedo model and Bennu's normal reflectance (Figure 16). The Lommel-Seeliger, Akimov, and Lunar-Lambert disk functions performed similarly to one another when paired with the ROLO and linear phase functions. Though the regional Akimov disk function had the lowest residual, it is not meaningfully better than the rest.
Even with a high-resolution shape model and regional disk function modeling, we were not able to accommodate the scatter in the data, as shown in Sections 3.2 and 3.3. Moreover, the variations that were apparent in the disk function maps (Figure 12) are likely influenced by uncorrected regional phase function variation. As such, for photometric correction, we prefer the Lommel-Seeliger disk function, which produces similar results without requiring a regional solution.
The residuals of the best models are similar; there is no statistical reason to change the conclusion from Golish et al. (2021a) that Bennu's surface is best approximated by a global Lommel-Seeliger disk function with a regional ROLO phase function. The linear phase function shows some improvement over ROLO, but its lack of an opposition surge term has positive and negative consequences. Without an opposition surge term, photometric correction to normal conditions (0°, 0°, 0°) will always be underestimated by ∼15%, as shown in Golish et al. (2021a). However, because we cannot model the ROLO Figure 16. Photometric albedo maps, at 0.55 μm, for each disk function solved with a linear phase function.

Phase Function Disk Function Residuals
L-S Minnaert (Global) Minnaert (Regional) L-L (Global) L-L (Regional) Akimov (Global) Akimov (Regional) opposition surge term regionally, the implementation of a global opposition surge may be incorrect. Though the available data do not make it possible to quantify to what extent the opposition surge might vary on Bennu, we expect that it is no more than a ±15% variation (the approximate variation of phase slope on Bennu) on a ∼15% effect. Nonetheless, we advise caution when choosing between these models. When correcting to standard laboratory conditions (30°, 0°, 30°), an opposition surge term is not necessary; we recommend use of the linear model. It is also useful to compare the regional solution directly with the previous global solution. We implemented the previous global solution in the regional photometry environment to evaluate the models with the same error metric. That is, we forced the solution to match the results published in Golish et al. (2021a), but we reevaluated the data using the regional analysis least-squares solver. Table 5 lists the residuals for each phase/disk function pair used in the previous global analysis (Golish et al. 2021a) and the regional analyses from this work, with and without inclusion of the  130 phase data (which the original global analysis did not include). We removed the Approach images from the global data set before making this comparison.
When including  130 phase data, the regional models perform 2-10× better than the global models. At very high phase, a regional analysis performs better, because the inclusion of the high-resolution shape model removes some disk function variation and masks shadows more effectively. When the data are limited to 90°phase, however, the global analysis performs better than the regional. We conclude that the averaging inherent to the global analysis (each "data point" was the median reflectance and photometric angles for an entire image) smoothed out much of the high spatial frequency variation on Bennu, at the expense of losing regional (photometric or albedo) variation. For phase angles under 90°, this trade seems favorable. To test this, we repeated the regional analysis using square bins with widths of  2 ,  4 , and  8 and found no substantive change in the residuals. Averaging larger swaths of Bennu's surface will combine the opposing effects of smoothing the disk function and ignoring true regional variations. Nonetheless, the regional analysis allowed us to include a wider phase angle range without loss of fidelity in the phase function or disk function models (Section 3.3).

Phase Slope Map
While regional modeling did not dramatically improve the model residuals, it does provide insightful photometric maps of Bennu's surface. Plotting the parameters for most models is not meaningful, as the parameters are interdependent. However, the linear phase curve, which is one of our preferred phase models, has a single parameter. Figure 17 is a map of the phase Table 5 Medians of the Squared Norm Residual for the Global Solutions Published in Golish et al. (2021a), at 0.55 μm, Compared to the Regional Solutions from This Analysis in Table 4 Phase/Disk Functions Global (0 < α < 95°)   slope-the linear model's β parameter when using a regional Akimov disk function (Tables 2; 3)-scaled relative to the average (∼−0.023). Though we recommended use of the Lommel-Seeliger disk function for photometric correction (Section 3.4), here we used a regional Akimov disk function to ensure that we did not misinterpret regional disk function variations as phase function variation. However, as mentioned in Section 3.3, the reverse is also possible-that we are misinterpreting phase function variation as disk function variation. Fortunately, although the choice of disk function changes the phase slope of the entire surface by 5%-10%, it has no substantive effect on the relative trends discussed below. This phase slope map represents the relative steepness of the phase function, i.e., how quickly the surface darkens with increasing phase angle. Phase slope is commonly used as a tracer for other surface characteristics, such as roughness, because increased roughness results in stronger shadows and a steeper phase slope (Schröder et al. 2017;Li et al. 2019).
We examined the spatial trends in the phase slope map by taking the median in the latitudinal and longitudinal directions, smoothing the result with a  5 kernel ( Figure 18). The spatial structure is very similar to the structure seen in the albedo map (Golish et al. 2021c). Regions of lower albedo, such as Roc Saxum (∼30°longitude) and Tlanuwa Regio (∼250°longitude) in the southern hemisphere, are also regions of steeper phase slope and drive the southern hemisphere as a whole toward steeper slopes.
When we overlay this structure directly on the normal albedo map (Golish et al. 2021c), we see correlation between the large dark boulders or boulder fields and regions of steeper slope (Figure 19(a)). Multiple studies have shown this same relationship between albedo, roughness, and phase slope for other surfaces (Buratti & Veverka 1985;Shepard & Campbell 1998;Shkuratov et al. 2005), though separating albedo from roughness is challenging. In Figure 19(b), we plot the phase curves (again using the linear phase function with a regional Akimov disk function) for several representative locations on Bennu's surface (marked in Figure 19(a)). Figure 19(c) plots the same data normalized by Bennu's average to better visualize their separation at high phase. The locations chosen include a range of surface types, from large, rough, dark boulders (Tlanuwa Regio, β=−0.0258 ± 0.004), to the largest boulder on Bennu that has a phase slope similar to its average (Benben Saxum, β=−0.0228 ± 0.004), to an expanse of relatively bright and smooth terrain around Bralgah Crater (Bralgah flow field; Perry et al. submitted; β = −0.0199 ± 0.004).
To explore these linked characteristics, we plot the albedo, as measured by Golish et al. (2021c) and downsampled to 1°×1°bins, against the phase slope ( Figure 20). Though the data are noisy (R 2 =0.425), there is a linear relationship between the two. However, dark boulders on Bennu tend to have a rougher surface, whereas brighter boulders tend to be smoother and more angular (DellaGiustina et al. 2020). Therefore, this relationship does not necessarily break the albedo/roughness ambiguity.
Alternatively, if we examine metrics for roughness on the surface, we continue to see spatial correlation with the phase slope map. Figure  Here we smoothed the phase slope map with a 5°×5°median filter to reduce its resolution to match the spatial resolution of the OTES and OVIRS instruments (5°-8°spot sizes), and we plot the correlation between phase slope and thermal roughness in Figure 21(c). Although thermal roughness notionally measures roughness at larger spatial scales (∼1 cm) than photometric roughness, roughness variations tend to propagate down to smaller spatial scales owing to the fractal nature of planetary surfaces (Shepard & Campbell 1998). However, the linear correlation (Figure 21(c)) between phase slope and thermal roughness is weaker than between phase slope and albedo.
Similarly, Figure 22(a) overlays contours representing the standard deviation of the laser altimeter data used to derive the OLA v20 shape model. Here we use standard deviation of the data as a measure of macroscopic roughness at the scale of the OLA spot size (∼5 cm). We median-filtered the OLA data with a 1°×1°pixel kernel to match the spatial resolution of the phase slope map. Figure 22(b) shows the OLA data as the map overlaid with contours representing the phase slope map. As with thermal roughness, we see a spatial correlation between the photometric and macroscopic roughness, though the statistical correlation (Figure 22(c)) is also weak. OLA measures roughness an order of magnitude higher (~10 cm) than the thermal roughness, putting it further still from photometric effects. That the correlation here is weak is perhaps not surprising. Nevertheless, there are some noteworthy spatial associations. The many large (1-10 m) boulders in Tlanuwa Regio produce numerous sharp edges (rapid elevation changes) over short distances that result in high standard deviations in the OLA map. Meanwhile, the Bralgah flow field is a smooth region with fewer discrete large boulders and a correspondingly low standard deviation. On the other hand, Roc Saxum, which appears rough in the phase slope and thermal roughness maps (and visually), is only slightly rougher than the Bralgah flow field. This is perhaps because the Roc is a single large outcropping with a rough surface, but without many boulder edges.
We see more evidence of this relationship in an equatorial crater located at approximately  150 longitude (Figure 23).  There the rim and walls of the crater are obvious as regions of relatively low phase slope (Figure 17), whereas the center of the crater has a slightly steeper slope than Bennu's average. The slope in the center of this crater is very close to that of the crater in which Nightingale, the site of the OSIRIS-REx sample collection, is located. The annulus outside the crater's rim is populated by boulders and has a steeper median phase slope. Crater centers on Bennu are typically darker and redder (DellaGiustina et al. 2020), with finer materials, than average Bennu (Burke et al. 2021). That we find crater centers to have steeper slopes suggests either that the slope is driven by the low albedo or that the spatial scale of the photometric roughness includes the fine particles therein, both of which are consistent with our expectations. However, overlaying the median albedo and phase slope of the three regions of the crater on a global correlation plot (Figure 23(d)) demonstrates only a slight  albedo increase in the center of this equatorial crater. Other regions of very low phase slope along the equator uniformly correspond with walls and rims of other equatorial craters.
We have not been able to separate albedo and roughness with the available lines of evidence. We can conclude that phase slope is related to both albedo and roughness, as expected, but cannot determine the contribution from each factor. Photometric analysis of individual boulders and other surface features may help resolve this degeneracy if those features provide examples of dark and smooth or bright and rough surfaces.

Phase Reddening Map
Phase reddening manifests as spectral reddening with increasing phase angle (Gradie & Veverka 1986) and is commonly observed on airless bodies (Li et al. 2015), including Bennu (DellaGiustina et al. Golish et al. 2021a). We calculated the phase reddening by fitting a line to the phase slope terms as a function of the wavelengths of MapCam's color bands (Rizk et al. 2018;Golish et al. 2020). The map of phase reddening (Figure 24), however, is largely featureless, showing a gentle phase reddening across Bennu. The median phase reddening (calculated by remapping to an equal-area sinusoidal projection and taking its median) from 0.47 to 0.86 μm is 0.0013±0.009 μm −1 deg −1 . The 1σ uncertainty is high, in line with the noisy map, but the median is consistent with the previous global analysis (0.0017 μm −1 deg −1 ; Golish et al. 2021a). Moreover, it is consistent with several analyses using spectral data acquired by OVIRS over the same wavelength range (0.0014 μm −1 deg −1 ; Fornasier et al. 2020;Zou et al. 2021;Li et al. 2021), which have higher spectral resolution. It is slightly less than the phase reddening measured on Ryugu (0.002 ± 0.0007 mm −1 deg −1 ; Tatsumi et al. 2020) and much less than found on Ceres (0.0046 μm −1 deg −1 ; Ciarniello et al. 2017) and 67P (0.0104 ± 0.0003 μm −1 deg −1 ; Fornasier et al. 2015).
The phase reddening map shows longitudinal arcs that are artifacts, likely due to the primarily nadir viewing conditions at equatorial latitudes that result in high incidence and emission angles toward the limb. These artifacts were not visible in any other parameter map we explored. However, the phase reddening map is a synthesis of the entire data set, using all four color bands, and the differences between the filters are small. The phase reddening value for a given latitude/longitude bin is a fit to four data points. Relatively small observational biases were amplified throughout the processing pipeline, producing artifacts. Nonetheless, the lack of meaningful spatial structure over MapCam's wavelength range is consistent with a similar analysis using OVIRS data .

Impact of Roughness on Photometric Modeling
Our empirical modeling took into account a high-resolution shape model, phase angle variation across an image, and albedo variation across the surface, with regional disk and phase functions. Yet the empirical models only fit the surface successfully by averaging much of that surface together to reduce scatter. The regional phase functions are valuable to map albedo and phase curve parameters but are statistically improved over previous global models only when including data with very high phase angles (up to ∼130°).
These models have a proven history in representing a variety of surfaces, though none as rough as Bennu. Section 3.2 showed Figure 23. Zoom-in on the global basemap  shows an equatorial crater with relatively smooth walls and a rougher center (panel (a)), but the crater is not obvious in the normal albedo map (panel (b)). Boulders populate the area just outside the crater rim. These roughness variations are apparent in the phase slope map (panel (c)). A correlation plot between the albedo and phase slope for the three regions shows a slight albedo dependence (panel (d)).
such large variation in the data that we suspect the cause to be related to the data themselves. Our photometric data depend directly on the accuracy of the shape model and the registration of the images to that model. Both of these have been shown to be accurate at or below the scale of a pixel (25-33 cm pixel -1 ) in the images considered here  and Section 2.2). Therefore, we consider variations below the pixel scale. The highest-resolution images that OCAMS has acquired, during the Recon C campaign of the mission, had a pixel scale of~0.33 cm pixel -1 , approximately 100×higher resolution than the MapCam data used in this analysis. Even in the Recon C data, the surface continues to show roughness and shadowing at the subcentimeter scale, emphasized in a 250×250 pixel excerpt from an image acquired at 2020 March 03 20:57:00 UT (Figure 25). When we represent a surface such as Bennu's with 20 cm shape model tiles (red triangles) and image it with 33 cm pixels (orange square), we do not capture the full photometric variation. Subpixel shadows (microshadows) and emission/incidence angle variation likely affect the resulting reflectance of that pixel in a way that cannot be modeled without higher-resolution photometric data.

Hapke Modeling
Hapke's photometric formalism (Hapke 1981(Hapke , 1984(Hapke , 1986(Hapke , 2012b was developed, in part, to model planetary surfaces based on physical principles, rather than empirical fits. In particular, Hapke's roughness formulation (Hapke 1984) addresses subpixel variation of incidence and emission angle, leading to the development of "effective" incidence and emission angles and a roughness term that can affect the reflectance of a surface. This is well suited for Bennu's surface given the subpixel roughness (i.e., incidence and emission angle variation) illustrated in Figure 25 Here we use notation from Hapke (2012a); μ 0e and μ e are the cosines of the effective incidence and emission angles, respectively, corrected for roughness. These are intrinsically linked to the derivation of the surface roughness correction, q a q S , i, e, , where (¯)¯is the mean slope angle of the unresolved facets imaged by a pixel.
Two parameters in this formulation are interdependent: the single-scattering albedo (SSA), w, and the single-particle phase function (SPPF), p(α). The latter describes the scattering of light by individual particles as a function of phase angle. As we will see in our uncertainty analysis, these two terms trade against one another. To attempt to separate these factors, we explore the SPPF in two forms: single-parameter Henyey-Greenstein (1pHG; Equation (9)) and two-parameter Henyey-Greenstein (2pHG; Equation (10)). The 1pHG form is dependent on the asymmetry factor, ξ, which varies from completely forward scattering (ξ = 1), through isotropic scattering (ξ = 0), to completely backscattering (ξ = −1). The 2pHG form is similar but breaks the asymmetry factor into  two terms (ξ = −bc), which necessarily have valid ranges of 0 <b< 1 and −1 <c< 1: x a x x a x = -+ * * + p , 1 1 2 cos 9 2 2 3 2 The Chandrasekhar H-function, H(μ, w), approximates the multiple scattering that results from isotropic scattering and takes many forms in the literature. We adopt the following form, described in Hapke (2002): , and x is either m e 0 or m e . B SH represents the shadow-hiding opposition effect (SHOE), defined by its amplitude (B 0 ) and its width (h): As mentioned in Section 2.1, Bennu subtended a few dozen pixels in our only low-phase-angle image data (from the Approach mission campaign), preventing us from modeling the opposition effect regionally. As such, we adopted the SHOE parameters derived from the disk-integrated analysis in Hergenrother et al. (2019), B 0 = 2.06 and h = 0.11. This approach has the same weakness as the ROLO phase function analysis, as we had no ability to model the opposition surge regionally. We did not include a coherent backscatter opposition effect (CBOE) term in this analysis, as CBOE is typically associated with multiple scattering, which is expected to be minimal on a low-albedo surface.

Global Hapke Analysis
Because the previous global analysis did not include a Hapke solution, we first fit this model to the global photometric data described in Golish et al. (2021a), which treated Bennu as a uniform surface and only included images with phase angles <95°. This Hapke formulation has either three or four free parameters, depending on the choice of SPPF, because the SHOE parameters are held fixed. Nonetheless, the parameter space is nonlinear and interdependent. We used a constrained, derivative-free MATLAB solver, fmindsearchbnd (D'Errico 2020), to minimize the error of the solution, which we measure as the rms of the difference between the measured and modeled RADF. This algorithm uses the Nelder-Mead simplex search method (Lagarias et al. 1998), which iteratively expands and contracts its search region until it finds a minimum. Table 6 lists the Hapke parameters for this global analysis. The 2pHG version predicted somewhat higher SSAs (and, correspondingly, a lower asymmetry factor, ξ = -bc) than the 1pHG model. We calculated the geometric albedo for each filter to produce a Bennu spectrum ( Figure 26). Both models generally agree with the ground-based Bennu spectrum (Clark et al. 2011). The ground-based spectrum has been normalized to match the normal albedo as measured by OCAMS, which depends on its radiometric calibration (Golish et al. 2020), enabling a valid filter-to-filter comparison.
The rms errors for the 1pHG and 2pHG versions are similar enough that neither stands out as superior. Similarly, plotting the modeled RADF values as a function of the measured RADF shows that both adhere to an identity line ( Figure 27). Plotting the ratio of the two (measured/modeled) as a function of phase, incidence, and emission angle is a better measure of where the model fails to predict the surface behavior ( Figure 28). At high photometric angles, the model has as much as ∼35% errorsimilar to the amount of scatter seen in the empirical models (Section 3.2).
The search algorithm did not guarantee a global minimum but operated reliably if seeded with reasonable starting values from a coarse search. Figure 29 depicts the 1pHG Hapke parameter space for the global photometric solution, using MapCam's v filter. As mentioned above, the SSA and asymmetry factor are interdependent. Though there is an optimal solution, there is a relatively broad range over which shifting w can be accommodated by matching ξ. While this is expected behavior, it is less expected that the model has an extremely weak dependence on q. Despite, or perhaps because of, Bennu's extremely rough surface, the Hapke formulation slowly converges on a specific roughness parameter.
To further explore the sensitivity of these parameter solutions, we performed an uncertainty analysis modeled after Li et al. (2019). We fixed one of the Hapke parameters at a given value and solved for the other parameters using the techniques described above. We repeated this for a range of values for the fixed parameter and calculated the rms error at each point, producing a sensitivity plot for each ( Figure 30). We plot the optimized value of the fixed parameter and the values at which the rms error increases by 50%. The sensitivity analysis shows that the SSA and asymmetry factor are somewhat constrained over a physically realistic range of values. The roughness parameter, however, can vary considerably (±10°). The 2pHG model is even less constrained ( Figure 31). The allowed variation (again for a 50% increase in rms) in the SSA is much larger, and the roughness parameter is virtually unconstrained between  15 and  30 . The b parameter has a symmetric variation to the SSA, also broad. The c parameter is constrained to positive values, indicative of backscattering, and agrees with the negative ξ values in the 1pHG solution. However, among positive values, c is essentially unconstrained.
These sensitivity analyses suggest that the Hapke model parameters are not well constrained by our photometric data. However, because the 1pHG is the more sensitive of the two, and because their global residuals were roughly equal, we elected to use the 1pHG model for our regional Hapke analysis.

Regional Hapke Analysis and Mapping
We applied the same techniques described in Section 4.1 to perform a regional Hapke analysis of Bennu's surface. We divided the surface into 1°×1°bins and fit each bin to our Hapke form using the same constrained search method. As expected, the amount of data in a given latitude/longitude bin was much less than that used in the global Hapke analysis, such that the relative noise within a bin increased. Figure 32 depicts the maps for the SSA, asymmetry factor, and roughness parameter. The sensitivity maps (Figure 33) depict the full width of the sensitivity, calculated as shown in Section 4.1, which can be asymmetric. The median sensitivities for SSA, asymmetry factor, and roughness parameter were +0.067/-0.010, +0.24/-0.064, and +28°/-20°, respectively. The spatial structure in the map of rms error ( Figure 34) correlates with the structure in the sensitivity plots. Regions of lower rms error (higher sensitivity) for w and ξ are along the rims and walls of equatorial craters and in the smooth flow field north of Bralgah Crater (Figure 19). The flow field has the highest roughness sensitivity (lowest sensitivity width) on the surface.
This sensitivity should not be treated as a Gaussian distribution with those widths, nor should we conclude that the weakly constrained parameters are meaningless. The geometric and Bond albedo maps ( Figure 34) describe Bennu's albedo distribution well. The R 2 value of the correlation is very high (0.986; Figure 35(a)). Moreover, the SSA and asymmetry factor have similar maps that follow the distribution of Bennu's surface features, including those highlighted in Figure 19. A correlation plot (Figure 35(b)) depicts the expected trade-off between SSA and asymmetry. However, it also shows that the  analysis predicts a broader range of SSA and asymmetry than is physically reasonable. Though Bennu's surface does have a small population of bright (albedo >0.07) features DellaGiustina et al. 2021), they are not nearly the ∼4% of the surface predicted by this analysis. These unphysical values correspond to the high-frequency noise (∼red pixels) in the maps, which are most prevalent in Tlanuwa Regio (dominated by large, dark boulders) and areas with a high density of medium-sized bright boulders . These regions are observed with generally high emission and incidence angles owing to the density of large boulders. For the same reason, they are the regions that have high standard deviations in the OLA roughness map (Figure 22). It is likely that these increase the noise in their Figure 28. Plotting the measured/model RADF ratio, at 0.55 μm, as a function of phase, incidence, and emission angle shows significant scatter at higher photometric angles. Figure 29. Global Hapke analysis of Bennu at 0.55 μm shows a parameter space that is weakly dependent on q and interdependent between w and ξ.
latitude/longitude bins, resulting in aberrantly high SSA values. These regions also correspond to areas of the sensitivity maps with higher uncertainty.
Any structure that might exist in the roughness parameter map ( Figure 32) is less obvious. There is an apparent equatorial band, though it is subtle-having a median value only ∼5°higher than the midlatitudes. The parameter sensitivity maps ( Figure 33) show that the roughness parameter was never constrained by less than ±5°, and more often by much more. There are also vertical bands of high roughness along the equator that are reminiscent of the longitudinal arcs visible in the phase reddening map (Figure 24). These are clearly observational artifacts and should not be interpreted as physical. High roughness bins (red speckles in Figure 32) other than the vertical bands are typically associated with high densities of boulders, i.e., the regions that have high SSA values. Nonetheless, the roughness parameter map appears to have a reasonable signal-to-noise ratio; the Hapke analysis consistently fits most of the surface to a range between  15 and  30 . Moreover, other analyses of Bennu have identified distinctive characteristics of Bennu's equator, including topographical , geologic , thermal (Rozitis et al. 2020), and spectral Simon et al. 2020;Li et al. 2021) properties. The equator is a geopotential low on Bennu (Scheeres et al. 2019;Barnouin et al. 2020), causing it to accumulate material moving downslope from midlatitudes ). As such, midlatitudes have large, exposed boulders, whereas the dearth of such boulders at the equator suggests that they may have been buried under deposited material ). This is seemingly at odds with the higher photometric roughness at the equator indicated by our analysis, although the photometric roughness scale is orders of magnitude smaller than the geologic roughness. Figure 36(a) plots the median roughness (smoothed with a  5 kernel) as a function of latitude; the equatorial difference is pronounced. The longitudinal cross sections for the northern and southern hemispheres (Figure 36(b)) indicate that they are largely similar, except for regions of high boulder density.
Bennu, at least on a 1°scale, may simply not have significant regional variation in its photometric roughness, as represented by the Hapke model. There is no correlation between the albedo map (Golish et al. 2021c) and the Hapke roughness parameter (Figure 35(c)). In Section 3.5, we noted that both albedo and thermal roughness (spatial scale of ∼1 cm) correlate with phase slope. Yet our Hapke analysis, theoretically measuring roughness at the photometric scale (∼100-1000 μm; Helfenstein & Shepard 1999;Hapke 2012a), seems to indicate no correlation with the phase slope. At face value, this suggests that SSA, rather than roughness, may impose a stronger effect on phase slope variations on Bennu. However, because the sensitivity of the Hapke roughness map is low (Figure 33), it may instead tell us that the Hapke model, as implemented here, does not represent the photometric roughness well. In that case, it is interesting that the regions in the roughness map that are more sensitive (lower values)  include the Bralgah flow field, the area southwest of Tlanuwa Regio, the area southwest of Roc Saxum, and the rims/walls of equatorial craters. These relatively smooth regions are less noisy and might therefore be easier to fit to the Hapke model. All of these regions have shallower phase curves (Figure 17). Whereas the absolute values of the Hapke roughness map may lack structure, the sensitivity plot may indirectly confirm that relatively smooth regions, which are easier to model, also have relatively shallow phase slopes.

Conclusions
Our regional photometric analyses have demonstrated Bennu to be an extremely challenging surface to model. Neither a suite of empirical models nor a Hapke model were able to represent the variations across the disk faithfully, particularly at high photometric angles. We conclude that the solution recommended by Golish et al. (2021a)-a ROLO phase function paired with a Lommel-Seeliger disk function-performs as well as or better than the other empirical models tested in this Figure 32. Hapke 1pHG parameter maps at 0.55 μm reveal the expected structure of Bennu in SSA and asymmetry factor, which correlate with albedo. The roughness map is largely featureless; the apparent equatorial band is weakly constrained and contains artifacts. regional analysis. However, we were unable to model the ROLO opposition surge owing to lack of high-resolution, lowphase data, and therefore we use the global opposition surge parameters. When correcting to standard laboratory conditions (30°, 0°, 30°), where an opposition surge correction is unnecessary, we recommend use of the linear model.
The difficulty of the photometric modeling is due to Bennu's high degree of roughness, which is more dramatic than most studied planetary surfaces, with the exception of Ryugu (Sugita et al. 2019). In the global analysis of Ryugu presented in Tatsumi et al. (2020), ∼60% of the data fit within 5% of the Hapke models (that is, modeled I/F values were within 5% of the measured I/F values). However, Tatsumi et al. (2020) used only phase angles (α) up to 50°. By analyzing Bennu's surface regionally, and with a high-resolution shape model, we find that ∼90%-95% of the α<50°data in our Hapke analysis fit within 5% of the model. Approximately 60%-65% and 50%-55% of the data are within 5% of the model when we include phase angles up to 90°and  135 , respectively-a dependence that is expected and exaggerated on an extremely rough surface. Nonetheless, the high spatial resolution of the image data and the unprecedented resolution of the shape model permitted modeling of even the high-phase data.
In addition, these regional analyses deepen our understanding of Bennu. The variation of phase slope of the linear phase function (β) across Bennu's surface correlates with albedo and roughness, which in turn correlate with each other. This is in line with our expectation for most planetary surfaces, but it is interesting on Bennu owing to its high surface roughness. We find that albedo is a slightly stronger indicator of phase slope variation than independent measures of roughness, but that both correlate statistically. However, the analysis of an individual equatorial crater showed a distinct change in phase slope without an obvious change in albedo. Unfortunately, the lack of structure in the Hapke roughness map does not help to confirm this, presumably because it is weakly constrained. Nonetheless, the relative sensitivity of the Hapke roughness does correlate with shallower phase slopes, Figure 34. The Hapke 1pHG geometric albedo and bond albedo maps at 0.55 μm match the published Bennu albedo map (Golish et al. 2021c). The rms error map, also at 0.55 μm, shows spatial structure similar to the sensitivity plots ( Figure 33). higher albedo, and less macroscopic roughness. Moreover, both the Akimov and Lunar-Lambert regional disk function analyses indicate that macroscopically rough regions correlate with high photometric roughness. Photometric analyses of the Moon (Kreslavsky et al. 2000;Velikodsky et al. 2000Velikodsky et al. , 2002 showed that high Akimov disk function roughness was directly proportional to surface roughness and inversely proportional to albedo. These multiple lines of evidence do not converge on a single answer, but they suggest that both albedo and roughness are correlated with phase slope. The dependence of phase slope on wavelength, over MapCam's four color bands, measures the magnitude of Bennu's phase reddening (∼0.0013 μm −1 deg −1 ) and is consistent with previous results Zou et al. 2021;Li et al. 2021). There is no coherent spatial structure to the phase reddening map, indicating that it occurs uniformly over the surface. Variations in its magnitude may be smaller than the noise in our measurements.
Sections 3.2 and 3.3 showed that our photometric data are very noisy, even after registration to a 20 cm GSD shape model, and in Section 3.7 we showed that the surface is still extremely variable beneath the 20 cm facets. Consequently, we conclude that the data, specifically the photometric angles, are not high enough in resolution to represent the reflectance variations measured by MapCam. In addition, the models used here may not be sufficient to describe Bennu's rough surface, particularly at high photometric angles. Hasselmann et al. (2020) applied a seminumerical model (van Ginneken et al. 1998) to describe the roughness component of the bidirectional reflectance distribution from MapCam Equatorial Stations data and found that a specular contribution improved the scattering behavior for emission angles higher than  50 . The origin of this specular contribution was not resolved in the images, but because the specular behavior was not represented in the models considered here, it may have contributed to our mismatch between model and data. Perhaps more importantly, DellaGiustina et al. (2020) showed that Bennu's surface includes at least two, and likely more, distinct types of boulders, even without including the very sparse, ultrabright exogenous surface material ). These boulder types on the surface are extremely well mixed (DellaGiustina et al. 2020;Golish et al. 2021c) at all scales from global to below the pixel resolution of the photometric data set. The regional analysis presented here does not have sufficiently high spatial resolution to separate the two main populations of boulders identified in DellaGiustina et al. (2020). However, future studies of individual boulders and boulder types (rather than binning the surface into arbitrary latitude and longitude bins) might reveal photometric differences between them and better discriminate regional characteristics of Bennu's surface.  (Golish et al. 2021c) is strongly correlated with the geometric albedo calculated in our Hapke analysis (panel (a)). SSA is also strongly correlated with asymmetry factor (panel (b)). Plotting the normal albedo against the roughness parameter shows no correlation (panel (c)).