Comparative photometric analysis of the Reiner Gamma swirl and Chang’e 5 landing site

Context. Lunar swirls are bright albedo features only found on the Moon that are still not entirely understood. It is commonly accepted that reduced space weathering plays a role in explaining the origins of lunar swirls because the local magnetic fields that are typically associated with these albedo anomalies are effective in reducing the solar wind influx. However, additional processes are required to fully explain the spectral, photometric, and polarimetric properties of the swirls. Aims. In this study, we compare the photometric properties of the Chang’e-5 landing site to those of the Reiner Gamma swirl. Because the physical effects of a landing rocket jet on the lunar regolith are relatively well known, these observations can provide important insights into the physical properties of lunar swirls. Methods. We determined the single scattering albedo, opposition effect strength, and surface roughness of the Reiner Gamma swirl and the Chang’e-5 landing site with their respective statistical uncertainties based on the Hapke model and Bayesian inference sampling. Results. The Chang’e-5 landing site and the Reiner Gamma swirl exhibit similar photometric properties, in particular: an increased albedo and a reduced opposition effect strength. Additionally, the landing site is about 20% less rough compared to the surrounding area. Conclusions. These findings suggest that the swirl surface is less porous compared to the surrounding surface, similarly to a landing site where the top layer of the regolith has been blown away effectively so that the compactness was increased. We conclude that external mechanisms that are able to compress the uppermost regolith layer are involved in lunar swirl formation, such as interactions with the gaseous hull of a passing comet.


Introduction
One of the biggest unsolved mysteries related to our nearest neighbor are lunar swirls. These swirls are characterized by a high-albedo pattern of complex shape and almost all swirls are located at a magnetic anomaly (e.g., Hood & Schubert 1980;Hood & Williams 1989;Tsunakawa et al. 2015). Furthermore, swirls have been found to show a reduced surficial OH/H 2 O content (Kramer et al. 2011;Hendrix et al. 2016;Li & Garrick-Bethell 2019;Hess et al. 2020a) and the spectral properties appear to be less mature than the surrounding terrain (Glotch et al. 2015;Hendrix et al. 2016;Trang & Lucey 2019;Blewett et al. 2021;Hess et al. 2020a). Previously, it was assumed that there is no correlation between topography and the albedo patterns (Kramer et al. 2011;Denevi et al. 2016), but more recent works have suggested that swirls are marginally lower than the surrounding terrain (Domingue et al. 2022). The photometry of swirls has also been observed to be anomalous. Kreslavsky & Shkuratov (2003) and Kaydash et al. (2009) showed that the phase function is less steep at Reiner Gamma and other swirls. Pinet et al. (2004) mapped the Hapke parameters with a genetic algorithm based on Clementine data and found a higher roughness and more forward scattering behavior for the Reiner Gamma swirl. Shkuratov et al. (2007) demonstrated that the Reiner Gamma swirl shows a weaker polarization than the surrounding mare. Bhatt et al. (2023) reported variations within the Reiner Gamma swirl structure from telescopic observations. They found larger grain sizes at the central oval compared to the extended tails but only marginal differences in surface roughness from its surroundings.
Several processes can contribute to the observed anomalies. One of the commonly accepted mechanisms is the shielding of the surface by the local magnetic fields, leading to a reduced solar wind flux (e.g., Wieser et al. 2010;Lue et al. 2011). The solar wind is one of the primary drivers for space weathering (Pieters & Noble 2016) as well as for the diurnal buildup of OH/H 2 O (Grumpe et al. 2019) in the uppermost layer of the regolith. However, the spectral trends at swirls are not entirely consistent with reduced space weathering. In the near-infrared (NIR) range, immature regolith is usually associated with a reduced spectral slope, but at swirls, this is not always the case (e.g., Pieters et al. 2021;Hess et al. 2020a;Blewett et al. 2021). In the ultraviolet (UV) wavelength range, the spectra are less blue (being spectrally blue refers to a relatively higher reflectance at a shorter wavelength) than would be expected for immature surfaces of similar composition (e.g., Hendrix et al. 2016). Therefore, additional processes are necessary to fully explain the occurrences of these swirls. It should also be noted that solar wind standoff does not imply a specific formation mechanism. To explain the unexpected spectral differences, a largely wavelength-independent brightening component is needed, where a difference in porosity is one of the possible explanations (Hess et al. 2020a). Garrick-Bethell et al. (2011) suggested that an electric field generated by the interaction of the magnetic field with the solar wind may sort small electrostatically charged dust particles. Furthermore, if swirls are enriched in feldspathic materials (e.g., Garrick-Bethell et al. 2011) or contain less TiO 2 , the higher albedo could also be explained. While Sato et al. (2017) found reduced TiO 2 content on swirls, the authors also noted that this is due to similar effects of maturity and TiO 2 abundances on the 321/415 nm reflectance ratio used for the estimation of TiO 2 . Another factor that might contribute to the optical properties of swirls is compaction (e.g., Hess et al. 2020a). Compact soil generally appears brighter but does not change the spectral slope and could possibly also explain the reduced phase curve steepness noted by photometric studies (Kreslavsky & Shkuratov 2003;Kaydash et al. 2009). Studies have suggested that the swirls have been (at least partly) formed as a result of impacts by a comet or meteoroids (e.g., Schultz & Srnka 1980;Hood & Williams 1989). The interaction of the surface with the comet's coma would also lead to the destruction of the highly porous fairy castle structure usually found on the Moon (e.g., Schultz & Srnka 1980;Shevchenko et al. 1994;Syal & Schultz 2015). Pinet et al. (2000) proposed a formation of swirls by impacts of fragments from a comet that has been torn apart by tidal forces. To be able to judge the relevance of the different processes, more information about the physical properties of the regolith at swirls is necessary.
Photometric analysis is an important tool to characterize the physical properties of planetary surfaces. However, imagery acquired under a multitude of different illumination and observation geometries is required (e.g., Shkuratov et al. 2011). Many studies are largely limited to the qualitative assessment of phase ratio images, which are difficult to interpret. To quantitatively evaluate the surface properties, knowledge about the topography and a suitable reflectance model are required. The parameters of semi-physical reflectance models, such as the Hapke (2012b) model, represent physical properties to some extent, but the parameters are oftentimes interrelated (e.g., Helfenstein 1986; so that effectively there is no unique solution. An entirely probabilistic perspective of the inversion problem, as in the case of Bayesian inference sampling, is well suited to account for these uncertainties and correlated parameters and such approaches have been successfully applied to, for instance, Mars (Fernando et al. 2013 or Jupiter's moons (Belgacem et al. 2020(Belgacem et al. , 2021. This decade can be seen as the renaissance of lunar exploration. There have not been as many missions aiming for the Moon since the Apollo era. In 2020, the Chinese lunar sample return mission Chang'e-5 successfully landed in northern Oceanus Procellarum, and the ascending stage lifted off again to return samples to Earth (Zhou et al. 2022). Landing sites of lunar missions have also been studied photometrically. Kreslavsky & Shkuratov (2003) investigated the Apollo 15 landing site and interpreted the less steep phase curve as a result of reduced roughness. Clegg-Watkins et al. (2016) observed similar effects for the Chang'e-3 and Apollo landing sites based on Narrow Angle Camera (NAC) images and concluded that the photometric properties are not significantly altered over the 40-50 years timespan between these missions. Xu et al. (2022) investigated the Chang'e-5 landing site at a very fine spatial scale with the lunar mineralogical spectrometer (LMS) and focused on the analysis of the phase function parameters that indicated a pronounced forward scattering behavior in the sampling zone of Chang'e-5.
The objective of this study is to compare the effects of a landing rocket jet, which are relatively well understood and correspond to a reduction of the roughness and destruction of the fairy castle structure of the regolith, with the photometric properties of the Reiner Gamma swirl for understanding the formation mechanisms of swirls. For the quantitative analysis and the investigation of the reliability of the data, we use Bayesian inference, which enables us to estimate the most likely solution along with precise values of the statistical uncertainties of the inferred model parameters.

Data sets and methods
For the photometric analysis of specific regions on the Moon, imaging from different illumination and observation configurations and accurate information about the geometrical layout of the scene are required. Therefore, the image data need to be processed with the Integrated Software for Imagers and Spectrometers (ISIS). Furthermore, suitable Digital Terrain Models (DTMs) are required to determine the surface orientation on the pixel level. In the following sections, we briefly describe the data sets (see Sect. 2.1) and the processing steps (see Sect. 2.2). After the geometry of each location is clearly defined, the Hapke model (Sect. 2.3) can be inverted using Bayesian inference to obtain the most likely solution and the associated uncertainties. We describe the Bayesian modeling in Sect. 2.4.

Lunar Reconnaissance Orbiter camera data
The Lunar Reconnaissance Orbiter (LRO) has been orbiting the Moon since 2009, carrying two cameras, which make it particularly interesting with respect to photometric analysis. The Narrow Angle Camera (NAC) returns high-resolution images of up to 0.5 m pixel −1 that can be used to identify boulders and can resolve rovers and landers on the lunar surface (Robinson et al. 2010). It has a narrow field of view (FOV) of 2.85 degrees, such that the phase angle is almost constant within a single image. A NAC image pair consists of a left and right image, which can be mosaicked using the 135 overlapping pixels. The sensor is sensitive to light in the visible wavelength range between 400 nm and 760 nm. The images used in this study for examining the surface before (see Table C.1) and after (see Table C.2) the landing of Chang'e-5 are listed in Appendix C. Before the landing, the area was imaged eight times with the LRO NAC instrument (Table C.1). Between the landing and the lift-off of the ascending unit, another image was taken that was excluded from this study. After the landing, another 11 images were taken of the landing stage, as listed in Table C.2. The images' phase angle range is between 46 and 79 degrees for the images taken before December 1st, 2020, and between 45 and 89 degrees after the landing. The region of interest was reduced to an area of 90 × 90 pixels or 144 × 144 m 2 size, in which the co-registration of the images was most accurate.
The LRO Wide Angle Camera (WAC) is a multispectral instrument that provides unprecedented coverage of the lunar surface from different phase angles. The camera covers a broad swath with 60 degrees FOV, such that the phase angle varies strongly within an image. The push-frame sensor A226, page 2 of 21 Hess,M.,et al.: A&A proofs, consists of seven wavelength channels between 321 nm and 689 nm (Robinson et al. 2010). Importantly, due to its design, each channel has a different emission angle. The WAC has a resolution of 100 m pixel −1 for the bands in the visible wavelength range and 400 m pixel −1 for the two channels in the near ultraviolet (NUV). In this study, only the 604 nm wavelength channel was used at a reduced resolution of 400 m pixel −1 . These WAC data have been studied extensively and various data products were derived, such maps of the Hapke parameters (Sato et al. 2014) binned into tiles of 1 × 1 degrees. We used the WAC images from the Experimental Data Record (EDR). The EDRWAC images of Reiner Gamma that we used are listed in Table C.3 for the oval part and in Table C.4 for the tail region.
To study the landing site of Chang'e-5, high-resolution image data were required. For this purpose, only NAC images provided the necessary resolution and sufficient number of images for different observation geometries. The Reiner Gamma swirl extends over a much larger area so that the WAC resolution was suitable and a wider range of phase angles could be accessed to better constrain the photometric model.
For the determination of photometric parameters, a highquality DTM is essential. Depending on the required resolution, a different generation method may be chosen. For the WAC images at 400 m pixel −1 , the GLD100 (Scholten et al. 2012) at 100 m pixel −1 was well suited. For the NAC images, no high-resolution DTM was available. Therefore, we used the Shape from Shading (SfS) framework described by Grumpe et al. (2014) to refine the photogrammetrically derived GLD100 to NAC pixel level resolution. The reference image for the co-registration of the NAC observations was M1114519536. Consequently, this image was also used to construct the DTM.

Data processing
To be able to use NAC and WAC images for estimating the photometric parameters on the pixel level, the images needed to be co-registered accurately, such that pixels with identical pixel coordinates in all images correspond to the same location on the ground. However, the images first needed to be radiometrically calibrated to radiances. For this purpose, we used the software ISIS provided by the US Geological Survey (USGS) (Laura et al. 2022). Due to the push-frame design of the WAC sensor, the images are split into even and odd images that need to be combined in a mosaicking step (the automos routine of ISIS). The images usually do not cover the entire area used for the map projection (Fig. 1a). The general steps required to process WAC images are displayed in Fig. 2. All even images were map-projected to the area displayed in Fig. 1 at a resolution of 400 m pixel −1 . The even image is then used for the map projection of the odd image. For the pixel grid of latitudes and longitudes, a list was created that was fed into the campt command in ISIS with the usecoordlist option set to True. Since each channel has a slightly different emission angle, the map-projected image cube was split up into the different wavelength channels with the ISIS explode command. The wavelength channel used in this study is the one at 604 nm. Therefore, the cube file for this wavelength channel was used for the campt command. The output for the even and the odd images was parsed to obtain the sub-solar and sub-spacecraft longitude and latitude values, as well as the solar distance and spacecraft altitude, for all the elements on the list. In cases where there was no sample taken for a specific image (even or odd), the values were set to "not a number" (NaN). The values were then reorganized in the grid (Fig. 1b) Fig. 2. Pipeline to radiometrically calibrate and geo-reference WAC images. The output is an image where each pixel location is defined in the coordinate system of the global WAC mosaic and the GLD100 DTM, and the angles i, e, and g are given for each pixel are based on the SPICE kernels. Purple denotes external data input. Orange denotes steps carried out in ISIS. Yellow denotes steps implemented in MATLAB or Python. was already almost perfect. Based on the latitude and longitude values of each pixel and the sub-spacecraft and sub-solar latitude and longitudes, as well as the solar distance and spacecraft altitude, which are also given for each pixel, we were able to calculate the vectors pointing to the sun (s) and to the spacecraft (u). The DTM is co-registered to the images, such that we could derive the surface normal n for each pixel. Subsequently, the vectors s, u, and n were used to calculate the angles i, e, and g. The resulting phase angle for an image showing a typical opposition surge is displayed in Fig. 1c.
The processing of NAC images ( Fig. 3) was easier regarding the ISIS calibration. The FOV is much smaller in comparison to WAC images. Therefore, the angles are not as highly variable as for WAC images. Also, no even and odd images needed to be combined. However, because of the higher resolution, much more attention had to be paid to the co-registration process. Consequently, all images were carefully co-registered to the reference image M1114519536, which was also used to create the SfS DTM. The image was much larger compared to the area of interest. Corresponding pixels for co-registration were selected, especially around the landing site, but also across the entire image. Furthermore, the alignment was evaluated and iterated until it was satisfactory. The SfS DTM is inherently co-registered to the reference image because the SfS method is directly based on the image information of that particular image. As with the WAC imagery, the vectors s and u were calculated based on the output of the campt command. Finally, the angles i, e, and g were obtained, which are required for the subsequent photometric modeling.

Hapke model
The Hapke (2012b) model is a widely used radiative transfer model for particulate planetary surfaces. In general, it is empirical, but the parameters are physically motivated such that they can be interpreted as physical properties of the regolith. For a detailed description of the model, see Appendix A, Hapke (2012b), and Sato et al. (2014). The parameters of the Hapke model are, on the one hand, either geometric (incidence angle, i, emission angle, e, or phase angle, g) or, on the other hand, physical (single scattering albedo, w, phase function parameters, b and c, opposition effect strength, B S 0 , and width, h s , as well as surface roughness,θ b ).
The opposition effect is related to the compaction of the regolith and leads to a strong peak for small phase angles; however, for typical lunar surfaces, this effect is still measurable at larger phase angles. The width of the opposition surge is modeled by the parameter h s , and the amplitude is denoted as B S 0 . A compact surface generally appears brighter compared to the highly porous fairy castle structure usually found on the Moon, but the opposition effect is less strong because there are fewer shadows that could disappear for small phase angles. This typical correlation between albedo and opposition effect amplitude is empirically described by Sato et al. (2014) for the Moon with  channel used in this work, s = 2.438 and t = 0.096. These values were also used for the NAC images because 604 nm is in the middle of the 400-760 nm sensitivity range of the camera. By dividing the estimated value by this predicted value, a corrected value can be obtained that gives insight into whether the estimated value is abnormal compared to the average lunar regolith, namely, smaller or larger compared to the average B S 0 for this albedo: (2)

Bayesian model
Certain parameters of the Hapke (2012b) model are interdependent (e.g., w and B S 0 ) and some are even mathematically coupled, such as the parameters of the phase function b and c (e.g., Helfenstein 1986;. Therefore, there are multiple solutions that can explain measured values equally well. The usage of conventional optimization techniques to invert the Hapke model may return best-fit values, but ignoring all other solutions that fit equally well. It is important to be aware of these uncertainties to be able to choose parameters that can be determined reliably based on the available image data. One technique that is well suited for these kinds of problems is Bayesian inference because the parameter space is thoroughly explored, and we obtain uncertainties complementary to the best-fit values. Because of the nonlinearities of the model and the not necessarily normally distributed parameters, we used Markov chain Monte Carlo (MCMC) sampling to numerically estimate the posterior probability density function (PDF) of the parameters, instead of solving the problem analytically. Choosing a suitable likelihood function is inarguably the most important part of the Bayesian modeling approach. The likelihood is evaluated by comparing the model results for a certain set of parameters to the measured data. If the modeled results and the measured data fit well, then it is likely that the data originated from these parameters. To be able to distinguish the effects of the albedo and the other parameters more easily, we additionally used phase ratio information. For M images, we defined the measurements as a vector, R, of the length of M. The first image is preferentially a medium phase angle image without any shadows. The phase ratios used in this work were the pixel-wise ratios R 1 /R 2...M between the reflectances of the first image divided by the reflectances of the remaining images and calculating the base-10 logarithm. This phase ratio vector has the dimension M − 1. Subsequently, we combined both by concatenation to obtain the final (2M − 1) × 1 data vector: Assuming Gaussian noise, it is reasonable to select a normal distribution as the likelihood function. The normal distribution we selected had its mean at the modeled values and was evaluated for the measured values. The standard deviation of the likelihood function is an additional parameter of the model itself.
Since the smaller the standard deviation, the higher the probability can become, and this has a strong effect. It should be noted that the reflectance values and the phase ratio values are of different orders of magnitude. Therefore, we included two different standard deviations as part of the model. We denoted the standard deviation for the reflectance part as σ refl and, for the phase ratio part, as σ pr . Thus, we defined the standard deviation vector σ lh of the likelihood function as: values individually and is subsequently marginalized. Using the definitions above, the likelihood becomes: In addition, we included information on priors. The Hapke model has many parameters and not all of them can be simultaneously determined in a reliable way. The single scattering albedo, w, which is the main parameter of the Hapke model, is physically limited to the interval between zero and one. We chose an uninformative prior with a Beta distribution evenly distributed between zero and one as: The opposition effects are collected in a correction term based on the width, h s , and the strength, B S 0 . The width is set to a fixed parameter according to Table 1. According to Sato et al. (2014), B S 0 is mainly distributed between 1.5 and 2.1. While in theory the amplitude is limited to the interval [0, 1] because both opposition effect terms are combined in B S , the amplitude is usually larger than 1. The parameter should not become negative so that we chose a log-normal distribution as a prior. In order to not overly penalize small values, we set the mean of the prior to e 0.4 ≈ 1.5 and the standard deviation to 0.7. Thus, we have: The surface roughness,θ b , has a very weak effect for roughness angles below 10 degrees and is rarely larger than 40 degrees. Therefore, we set the prior distribution to be equally distributed between 10 and 50 degrees by using a scaled Beta distribution: Besides the physical parameters of the model, the standard deviation of the likelihood distribution is also a free parameter. It is common practice to use a half-normal distribution for the prior of the standard deviation such that p(σ pr ) = HalfNormal(σ pr | σ hn = 1).
The parameters b and c of the phase function influence the shape of the phase curve and determine if the particles are forward or backward scattering and inherently follow the hockey stick relation (e.g., Hapke 2012a; , which is given by Sato et al. (2014) as c = 3.29e −17.4b 2 − 0.908. The parameter b is limited to the interval [0, 1]. Therefore, we again chose an uninformative Beta distribution as a prior for b.
To ultimately decide which parameters can be estimated, we designed an experiment where we estimated the uncertainties of the parameters given data synthetically generated with the Hapke model and a conservative estimate of the noise level (Gaussian noise with σ n = 0.001). The observational geometry was selected to represent the data available for the central oval of Reiner Gamma. Given the nine observations, each with the nominal angles i, e, and g, we calculated the reflectance values for w = 0.3,θ b = 23 degrees, B S 0 = 1.7, h s = 0.11, b = 0.21, and c = 0.619. The signal-to-noise ratio (S/N) is approximately 77.7 compared to >42 for NAC imagery, for example, not considering the radiometric accuracy. These reflectance values were then used to invert the Hapke model to determine the parameters by employing the Bayesian method outlined in this section. Figure 4 shows the results when the b parameter of the phase function is free and c is set according to the hockey stick relation (Hapke 2002). The results indicate that the phase function parameters cannot be determined reliably. All solutions between strongly backscattering and forward scattering are accepted. Comparing those results to the results when the phase function is fixed (see Fig. 5) visualizes that the uncertainties of all free parameters can be reduced significantly if the value for b is constant. Additionally, in situations where the mean deviates strongly from the actual value (yellow diamond) in Fig. 4 it fits well for Fig. 5. In practice, using b as a free parameter for the actual data often resulted in multimodal posterior distributions for w, b, and c, making the mean an unreliable measure to describe the distribution. As illustrated by the results, the parameter has only a limited influence on the reflectance values. The shape of the phase function simply does not change significantly enough to be able to distinguish variations in b. Figure 6 shows some examples of solutions that might possibly be accepted by the Bayesian inference sampling even for the very narrow confidence band. In this example, a strongly backwardscattering but low-albedo material fits the data equally well as a forward-scattering material with higher albedo because there is not enough angular coverage to constrain the model. The other parameters in the example in Fig. 6 are constant. This is an indication that variations in b and c between, for instance, the swirl and its surroundings, cannot explain the photometric differences we observe in this study. Therefore, we decided to fix the b and c values according to Warell (2004) and infer the values of w, B S 0 , andθ b from the data. The parameter choices and corresponding priors are listed in Table 1 posterior distribution. NUTS is a very efficient algorithm that is widely used in Bayesian applications. We used the implementation of NUTS in pymc3 (Salvatier et al. 2016) and sampled 5000 tuning steps to set the individual step size, drawing two chains of each 20 000 samples to make sure that the samples converge to the posterior density.

Results
In this study, we want to investigate two regions with respect to their photometric behavior and compare their characteristics.   (Zhou et al. 2022). This is an interesting study site because it offers the rare opportunity to investigate the photometric changes caused by a landing rocket jet because high-resolution NAC images from before and after the landing are available. Second, lunar swirls are still not fully understood and the Reiner Gamma swirl (see Sect. 3.2) is a prime research target as it is the most prominent representative of lunar swirls.

Chang'e 5 landing site
The Chinese lunar mission Chang'e-5 landed successfully on December 1st, 2020, at approximately 43.1 • N and 51.8 • W. It collected a number of samples and the ascending stage lifted off two days later and returned to Earth (Zhou et al. 2022). The context of the landing site is shown in Fig. 7, where the region of interest evaluated in this study is marked by the red rectangle. The area is characterized by a typical mare terrain with small craters scattered around. The landing site is salient by three small  Fig. 4 with some exemplary reflectance phase curves. Both strongly backscattering and forward-scattering material can fit into the confidence interval and could be considered acceptable solutions given the noise level. The shape of the phase curve only changes marginally across the accessible range of phase angles.
craters in a triangular pattern, where the craters are smaller than 10 m in diameter, as can be seen in Fig. 7b. The modeling results for image data acquired before the landing are shown in Fig. 8 as the mean predicted values (a-d), and the associated uncertainties are displayed as the standard deviation of the posterior distribution (e-h). The mean albedo of the area is 0.18, with no reliable values above 0.22. The roughness is overall small (mean over the area is 17.5 degrees) compared to the value in Sato et al. (2014), where a global mean value of 23.4 degrees was derived. In the maps of the roughness and the opposition effect amplitude, the craters are associated with higher values but are also associated with worse fits (Fig. 8d) and higher uncertainties (Figs. 8f and g). The mean B S 0 value over the area is 0.90 and, therefore, smaller than the values predicted for the derived low albedo. After the landing, the lander itself is visible in the albedo map as a bright spot between the three craters (Fig. 9a). The albedo increases significantly around the lander compared to before the landing (Fig. 8a). The surface roughness and the opposition effect are both reduced immediately around the lander, as can be seen in Figs. 9b and c in a similar pattern as the albedo has increased. Because of the Moon's orbit, the images are all illuminated from the east or the west but never from the north or south. Therefore, the shadows of the lander in the image are also cast in these directions. These shadows can be seen in the maps of the mean predicted parameters ( Fig. 9) as increased values in B S 0 andθ b and as higher uncertainties for all parameters as a bright line in the center of the image. Higher values of B S 0 andθ b , in general, are also associated with higher uncertainties, just as craters are (again) not as well constrained. The uncertainties in the craters also stem from shadows at the eastern and western edges of the crater, often leading to the two semi-circles with higher uncertainties and higher values visible in the B S 0 andθ b maps (see Figs. 8 and 9).
Using the correlation from Sato et al. (2014), the B S 0,corrected values are displayed in Fig. 10 for before and after the landing. A comparison of the two maps shows that the background has not changed significantly, but a strong negative anomaly can be observed around the lander.
Based on the albedo, which does not exceed 0.22 before the landing, we could segment the region of interest into an area that was affected by the landing (w > 0.22) and an area that was not (or at least less so). This segment acts as a background (w < 0.22). We also defined pixels for which the standard deviation of the likelihood function σ refl exceeds 0.004 as unreliable and omit them from our analysis. This segmentation is shown in Fig. 11.
Subsequently, this segmentation can be used to plot the histograms of B S 0,estimated and B S 0,corrected for the two regions (Fig. 12). In Fig. 12a, the B S 0,estimated values and in Fig. 12b the corrected values are displayed, both before the landing. The corrected values are generally below one (roughly normally distributed around approximately 0.53) even though we would not have expected this region to have an abnormally weak opposition effect. The histograms of the study sites show that the B S 0 values after the landing (Figs. 12c,d) are generally also small. The segmentation, however, clearly shows that the values of B S 0,corrected in the area affected by the lander, are excessively small, corresponding to less than one-third of the typical values of lunar regolith with a comparable albedo. The mean of the background pixels is also smaller, with 0.45 compared to 0.53 before the landing. The segmentation discretely classifies pixels as "landing site" or "background". However, the affected area actually has a diffuse boundary. Therefore, it can be stated that background pixels are also affected, even though to a lesser degree compared to the landing site. Consequently, the mean opposition effect strength of the background pixels is also slightly reduced compared to before the landing. The histograms of the two areas, landing site and background, appear to describe systematically different distributions. The mean of the B S 0,corrected values from the landing site is 0.28, compared to the background value of 0.45 after the landing. The Kolmogorov-Smirnow test (KS-test) statistics for the B S 0,estimated values of the two regions is D = 0.527 (p = 1.8 × 10 −15 ) and for the B S 0,corrected values it is D = 0.511 (p = 1.8 × 10 −15 ). The null hypothesis can thus be rejected.
The histograms for the roughness angleθ b are shown in Fig. 13 for before and after the landing. Before the landing, the roughness was overall higher with a mean of 17.5 degrees. After the landing, the background pixels appear to be only slightly affected by the landing rocket jet as the mean decreases to 17.2 degrees. The pixels classified as landing site show a reduced average roughness of 13.4 degrees. The modes of the two histograms, background and landing site, are very close to each other (12.5 and 11.6 degrees for background and landing site, respectively) and both are small compared to the mode before the landing, which was 14.7 degrees. The KS-test score is D = 0.436 (p = 1.78 × 10 −15 ), suggesting that the effect is not as large as for the opposition effect but still both regions are systematically different.

Reiner Gamma swirl
The Reiner Gamma swirl is centered at 7.5 • N and 59.0 • W in western Oceanus Procellarum, slightly north of the lunar equator. It is the most prominent example of these high-albedo anomalies found on the Moon that are commonly termed lunar swirls.
In the albedo map of the central oval part of the Reiner Gamma swirl, the bright albedo pattern is clearly visible (Fig. 14a). The roughness map (Fig. 14b) shows that the values are relatively noisy, and the uncertainties are high with several degrees and even exceed 10 degrees for some craters (Fig. 14f). A small difference can be observed for the dark lane in the A226, page 8 of 21 Hess,M.,et al.: A&A proofs,  northern oval, which exhibits a higher roughness. In the B S 0 map, however, the silhouette of the swirl can be distinguished by a weaker opposition effect amplitude. Craters, in turn, often coincide with weak B S 0 anomalies but are mostly related to higher uncertainties. The standard deviation of the likelihood function is higher for the craters (Fig. 14), indicating that the fit is worse, just as the uncertainties of all parameters are higher (Fig. 14). The results for the tail part are very similar to the central oval and are shown in the Appendix (Fig. D.1). The opposition effect strength is clearly smaller for the tail, even after correcting for the general correlation between albedo and B S 0 (Fig. D.2a). For the roughness in Fig. D.1b a very weak reduction can be observed at some locations, such as in the northernmost and southernmost parts of the tail.
When the dependence on albedo is removed from the opposition effect amplitude, the outlines of the swirl can still be clearly differentiated (Fig. 15a). The area can be segmented into onswirl pixels with a higher albedo and off-swirl locations that represent the surrounding background mare. The albedo threshold is set to 0.21 and locations with σ refl > 0.004 are labeled as invalid. With this segmentation, the histograms of the estimated and corrected values are shown in Fig. 16. The B S 0,estimated values for the off-swirl pixels are mostly normally distributed around 1.75, and for the on-swirl locations, the mean is 1.37. On the other hand, the off-swirl B S 0,corrected values are almost exactly centered around 1.0, which means that they are on average exactly as predicted by Sato et al. (2014). In contrast, the on-swirl locations have a reduced opposition effect with an average B S 0,corrected value of 0.83. Because of the large number of locations in the distributions, the KS-test is significant with p = 0.00 for both corrected and estimated values. For the estimated values, the statistics is D = 0.539 and for the corrected values the difference is smaller between the on-swirl and off-swirl distributions. Consequently, D = 0.422 is also smaller. This dichotomy between on-swirl and off-swirl locations is even more pronounced for the tail region of Reiner Gamma. The segmentation for the tail is shown in Fig. D.2b, and the resulting histograms are in Fig. D.3. For that part of the swirl, the KS-test statistics are D = 0.725 and D = 0.607 for the estimated and corrected values, respectively. The difference between the roughness values on-swirl versus off-swirl is negligible with D = 0.075 (p = 3.98 × 10 −33 ) for the central oval and D = 0.071 (p = 6.88 × 10 −9 ) for the tail.

Discussion
The bright regions of the Chang'e-5 landing site and the Reiner Gamma swirl show similar photometric behavior. Compared to the immediate surroundings, the albedo is increased, and the opposition effect is weaker in both cases. Additionally, at the landing site, the roughness is significantly reduced compared to the pre-landing surface, whereas a reduction of swirl roughness is barely discernible. The differences in opposition effect are similar for the landing site and swirl when compared to their respective surroundings, but the absolute values are excessively small at the landing site.
The strongly reduced B S 0 values at the landing site might be attributed to unsuitable standard parameters considered for the region. For example, the width, h s , of the opposition effect could be smaller at higher latitudes, but according to the maps of Sato et al. (2014), the width is almost constant within western Oceanus Procellarum (approx. 0.05 at both locations). Xu et al. (2022) found that the regolith at the landing site is more forward scattering compared to typical lunar regolith, which could also influence the estimated absolute B S 0 values. Since, for the landing site, only images at medium-to-large phase angles are available, the opposition effect is not fully constrained, as can be seen by the mean standard deviation of 0.7 over the landing site area before the landing and 0.51 after the landing, compared to 0.22 for the oval part of Reiner Gamma. The relative differences are still, however, strongly pronounced for the landing site. Importantly, the changes from the before to the after images show that a physical modification of the regolith structure has occurred.
This supports the findings of previous studies that the landing of a rocket on the lunar surface disturbs the microstructure of the regolith and leads to the destruction of the fairy castle structure, effectively making the regolith more compact (Kreslavsky & Shkuratov 2003;Clegg-Watkins et al. 2016). Xu et al. (2022) investigated a very small area around the Chang'e-5 lander and set the opposition effect and roughness parameters to constant A226, page 9 of 21    Fig. 11. Segmentation of the region of interest after the landing into an area strongly influenced by the rocket jet (landing site) and not as clearly influenced by the landing (background). Areas with a mean predicted albedo w above 0.22 are denoted as landing site. Areas with a mean predicted σ above 0.004 are labeled as invalid because the uncertainties are too high. Consequently, these areas are excluded from the analysis of the landing site and background regions.
values. Therefore, the results presented in this work are not directly comparable to Xu et al. (2020). Whether the observed behavior is due to a reduced shadow hiding opposition effect or more forward-scattering material cannot be reliably determined with NAC imagery. To effectively constrain the b and c parameters of the phase function, images acquired at phase angles well beyond 90 degrees are needed.
The presence of craters at both the landing site and the swirl coincide with higher uncertainties. The shape of the high uncertainty regions is not round as in the case of the crater but often split into two areas, one to the east and one to the west. The source of these uncertainties are shadows in large phase angle images. Considering the topography at each crater location alone, there should be light reflected even at the large phase angles. However, because of the crater topography, these regions are shadowed so that the model cannot adequately explain the large phase angle images, leading to higher uncertainties. In a qualitative analysis using only phase ratio images or in a modeling study with a classical least squares optimization, these areas would have to be carefully removed by hand, but because we are aware of the uncertainties, we can automatically omit them from the subsequent analysis.
The Reiner Gamma swirl consistently shows a reduced opposition effect strength, which is commonly interpreted to correspond to a more compact surface. At the landing site, we observed a similar behavior, and we know that the regolith has become more compact due to the rocket jet that blows away the porous top layer of the regolith. When comparing the Chang'e-5 landing site with the Reiner Gamma swirl, we find that both show a difference in compaction between the bright areas and the background. The difference found in opposition effect strength amounts to more than just the variation associated with albedo differences, suggesting that there must be a physical difference between the surface microstructure of bright and dark areas. The mean predicted roughness at the Reiner Gamma swirl shows only very weak anomalies correlated with the swirl albedo, whereas the landing site shows a more strongly reduced roughness. Shkuratov et al. (2010) theorized that the influence of multiple scattering of immature soils might counteract roughness effects. The gas pressure at the landing site is approximately 170 Pa on average in the exit plane and around 420 Pa at the nozzle wall (Zhang et al. 2022). At the swirls, Starukhina & Shkuratov (2004) calculated that a passing cometary coma creates pressures in the broad range of 40-400 Pa at the surface. Both regions exhibit pressures of a similar order of magnitude. The reduced roughness at the landing site is only visible in a small area around the lander of less than 20 m from the lander, where the pressure from the exhaust was highest. Possibly, the pressure of a cometary coma is closer to the lower end of the range and, therefore, not high enough to completely remove a regolith layer and transport it over much larger distances, but only to compact it without a significant change of the surface roughness. Additionally, the Chang'e-5 mission included an ascending stage that disturbed the surface structure of the regolith for a second time. Despite the presumably slightly higher exhaust pressure, the decrease in roughness is as small as 4 degrees (around 20%) or less at the landing site compared to before the landing.
Three mechanisms might be considered to play a role in the appearance of lunar swirls: (1) reduced or abnormal space weathering due to magnetic shielding (e.g., Trang & Lucey 2019;Blewett et al. 2021); (2) sorting of grains by electric or magnetic fields (e.g., Garrick-Bethell et al. 2011);and (3) interaction of a cometary coma with the surface (e.g., Schultz & Srnka 1980;Pinet et al. 2000;Starukhina & Shkuratov 2004  between the swirl and its surroundings. Therefore, a difference in compaction or grain size should still be considered a contributing factor. The individual contributions of the various processes might be different depending on the area. The spectral significance of compaction defined by Hess et al. (2020a) is higher in the northeastern tail region of the Reiner Gamma swirl, compared to the central oval. This is also consistent with the results presented in this work, where the difference between the opposition effect strength of the swirl and the surrounding terrain is more pronounced for the tail region.
A226, page 13 of 21 A&A 674, A226 (2023) N N N N N N N N N N N N N N   Also, the grain size strongly contributes to the observed spectral and polarimetric properties of the lunar regolith (e.g., Bhatt et al. 2023), but the grain size variation cannot be constrained based on photometric data alone. If the surface at swirls was enriched with finely grained feldspathic dust by electrostatic sorting as proposed by Garrick-Bethell et al. (2011), the spectral effect would be similar to higher compaction and could lead to a reduced roughness as well. However, the forces might not be strong enough to significantly alter the surface structure to have an influence on the surface roughness. While most swirls are colocated with magnetic anomalies, there is a swirl structure in the Mare Moscoviense basin that is not protected by a magnetic field (e.g., Hess et al. 2020a) that would be required to create the electric field responsible for the sorting of dust. A general swirl formation mechanism has to be able to explain these occurrences as well.
The photometric observations favor the interpretation that swirls are created by some interactions between a comet's coma and the lunar surface (Shevchenko 1993;Pinet et al. 2000;Syal & Schultz 2015). On the one hand, this could be attributed to the effects of the coma of an impacting comet. For example, Syal & Schultz (2015) showed that such a coma could cause scouring at scales of lunar swirls and can also explain the complex shapes. On the other hand, the passing of a comet and the interaction of the inner, denser parts of the coma with the surface would cause an effect similar to a landing rocket jet (Shevchenko 1993). The parts of a coma dense enough to be able to alter the physical surface properties are less than 1000 km in radius (Syal & Schultz 2015) and could explain the extent of Reiner Gamma, which has a length of roughly 300 km and a width of 45 km at the central oval. Magnetic shielding can preserve swirls by reducing the solar wind influx so that the surface does not become mature and more porous as quickly as the unshielded surface, which would mean that the Mare Moscoviense swirl not associated with a magnetic field is younger than other swirls (Hess et al. 2020b).

Conclusion
We introduce a Bayesian framework to estimate the Hapke photometric parameters and their uncertainties. Depending on the available data, a suitable model can be chosen and either of the following concerns can be directly inferred: which of the parameters are reliable, what additional data are required to improve A226, page 15 of 21 A&A 674, A226 (2023) the estimation accuracy, or which parameters are difficult to estimate or are mutually interrelated. In this way, the interpretation of the results can be focused on the areas with higher confidence.
We applied this framework to the Chang'e-5 landing site and the Reiner Gamma swirl. Both regions exhibit excessively weak opposition effects and an increased albedo. The reduction of the opposition effect amplitude goes well beyond the expected value for typical lunar regolith of similar albedo, suggesting that the physical properties of the regolith are also distinct. For the landing site (which additionally exhibits clear characteristics of reduced roughness), this optical change is a direct consequence of the top regolith layer having been blown away by the descent and ascent rocket jets, leaving behind a more compact and smoother surface. At swirls, this difference between the swirl and its surroundings could be attributed to changes caused by a cometary coma acting similarly to a landing rocket jet. There appears to be a minor difference between the roughness of the swirl and the surroundings, but the difference is statistically not significant when compared to the uncertainties. This behavior might be attributed to slight differences in gas pressure between the swirl and the landing site, where the pressure was presumed to be relatively high so close to the lander. A cometary coma might not generate enough pressure to alter the surface roughness, but it could still destroy the fairy castle structure of the regolith and increase its compaction. The magnetic shielding would then serve as a preservation mechanism resulting in immature surface regolith, with a fairy castle structure building up more slowly at the swirls than in the surrounding mare. . (A.10) The helper functions for both cases are defined as follows. The function χ is defined as (A.11) and the function η as The beta distribution is limited to the interval [0,1] and is defined as

B.2. Log-normal distribution
The log-normal distribution is given by: Lognormal (

B.3. Half-normal distribution
The half-normal distribution is a normal distribution with a mode at zero and all values below zero having zero probability.