Phase-curve analysis of comet 67P/Churyumov-Gerasimenko at small phase angles

Aims. The Rosetta-OSIRIS images acquired at small phase angles in three wavelengths during the ﬂy-by of the spacecraft on 9–10 April 2016 provided a unique opportunity to study the opposition effect on the surface of comet 67P/Churyumov-Gerasimenko (67P). Our goal is to study phase curves of the nucleus at small phase angles for a variety of surface structures to show the differences in their opposition effect and to determine which surface properties cause the differences. Methods. We used OSIRIS NAC images that cover the Ash-Khepry-Imhotep region to extract the phase curve, that is, the reﬂectance of the surface as a function of phase angle. We selected six regions of interest (ROIs) and derived the phase curves for each ROI. We ﬁt a linear-exponential function to the phase curves. The resulting model parameters were then interpreted by spectrophotometric, geomorphological, and phase-ratio analyses, and by investigating the inﬂuence of structural and textural properties of the surface. Results. We ﬁnd evidence for the opposition effect (deviation of the phase curve from linear behavior) in phase curves for all areas. We found an anticorrelation between the phase ratio and reﬂectance in a small phase angle range. This provides evidence for the shadow-hiding effect. We conclude that the decrease in the slope of the phase ratio versus reﬂectance indicates a decrease in the proportion of shadowed regions and reduces the contribution of the shadow-hiding effect. Large uncertainties in the determination of the opposition effect parameters with respect to wavelength do not allow us to conclusively claim coherent backscattering in the opposition effect phenomenon. Based on the two analyses, we conclude that the opposition effect of comet 67P in the Ash-Khepry-Imhotep region is mainly affected by shadow-hiding.


Introduction
The Rosetta spacecraft (Schulz et al. 2009) rendezvoused with its target comet 67P/Churyumov-Gerasimenko (hereafter 67P) in August 2014 and orbited the comet until 30 September 2016. During the 2.5-year mission, two zero-phase-angle fly-bys were performed by Rosetta. The phase angle α is defined as the angle at the comet between the Sun and the observer.
During the zero-phase-angle fly-bys, the scientific imaging system on board Rosetta, the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007), acquired high-resolution images of the comet surface in different filters in the visible wavelength range.
The first zero-phase-angle fly-by took place on 14 February 2015 with a closest-approach distance of 6 km from the nucleus surface. A study of this fly-by is presented in Feller et al. (2016) and Masoumzadeh et al. (2017).
The second zero-phase-angle fly-by took place on 9-10 April 2016. Rosetta reached a minimum distance of 30 km from the comet, and OSIRIS acquired images with the Wide Angle Camera (WAC) and the Narrow Angle Camera (NAC). The photometric analysis of the WAC images was performed by Hasselmann et al. (2017). The analysis in this paper is based on the OSIRIS NAC observations. At small phase angles, a phenomenon known as opposition effect (OE) manifests itself as a rapid increase in the surface brightness. The brightness dependence on the phase angle (known as phase curve) contains information about photometric and structural properties of surface. The OE can be described in terms of two parameters: the amplitude (also specified as the enhancement factor ζ) and its angular width, which is estimated as the half-width at half-maximum (HWHM).
Two mechanisms are believed to control the OE: shadowhiding (SH), and coherent backscattering (CB). The SH relates to the amount of shadows that particles cast on each other. The shadow cast by the regolith particles is progressively hidden from the observer as α approaches 0 • (Hapke 1986;Shkuratov 1994;Penttilä 2013). The second mechanism, CB, results from constructive interference between the partial electromagnetic waves that travel in the medium in opposite directions and experience multiple scattering on the same particles. These waves tend to leave the medium in phase and thus provide the conditions for the constructive interference. A&A 630, A11 (2019) Fig. 1. Footprint of 33 map-projected images (latitude range: −50 • to +30 • , longitude range: +70 • to +160 • ) superimposed on the cylindrical projection map of the shape model of 67P. The gray shading described the area on the shape model that was illuminated at 23:52:10 on 9 April 2016. The color gradient displays the overlapped area of the OSIRIS NAC images that were acquired during the fly-by on 9-10 April 2016. The color code is as follows for the overlapping: yellow shows all images, orange shows more than half of the images, bright blue is for less than half of the images, and dark blue shows none of the images. The map projection of the shape model was made with the ShapeViewer software (http://www.comet-toolbox.com/shapeViewer.html).
Our goal is to explore the photometric properties and microstructure of the surface of comet 67P, together with the physical mechanisms that can play a role in the OE phenomenon.

OSIRIS NAC data
For our study, we analyzed a dataset composed of NAC images acquired on 9-10 April 2016 that includes 99 images in three NAC filters. The NAC filters are F84 (central wavelength 480.7 nm), F82 (649.2 nm), and F88 (743.7 nm). The images were acquired when Rosetta was at a distance of 30 km from the nucleus center, corresponding to a spatial resolution of 0.53 m pix −1 . The images were acquired in triplets (one image per filter) within about 25 s at intervals of 3-6 min between triplets. As a result, the spacecraft (S/C) motion and rotation of the comet are small during the acquisition time of the triplet.
The NAC images cover an area in the Ash-Khepry-Imhotep region (El-Maarry et al. 2015) (see Fig. 1). During the fly-by, the phase angle to the center of the image decreased from 7.76 • to 0.99 • and then increased again up to 6.17 • . The observations started on 9 April at 23:01:09 UTC and ended on 10 April at 00:46:09 UTC.
In this paper we use OSIRIS level 3B images in radiance factor. The radiance factor, also known as I/F, is defined as the ratio between the bidirectional reflectance of an illuminated surface to that of a normally illuminated Lambert surface. We use the term reflectance for the radiance factor hereafter. The OSIRIS images were calibrated as discussed in the paper by Tubiana et al. (2015).
The absolute calibration has an uncertainty of 1-2% for NAC filters in the visible range.
For the purpose of photometric correction, we used the local scattering angles that were calculated with the global shape model of the comet (Preusker et al. 2017). The Lommel-Seelinger disk function (LS), which describes the photometric behavior of the surface of 67P sufficiently well, was used for the photometric correction (Fornasier et al. 2015). The photometric correction and coregistration of the images in three filters were made with the USGS ISIS software (integrated software for imagers and spectrometers 1 , Anderson et al. 2004). For our study, we selected six regions of interest (ROIs), shown in Fig 2. We chose these ROIs based on different geomorphological appearance (see Sect. 2.2) and because they were observed at very small phase angles (α < 1 • ), as shown in the center of Fig. 2.

Geomorphological analysis of ROIs
In this section, we classify our ROIs according to the geomorphological units on the Ash-Khepry-Imhotep region. The ROIs morphologically comprise three classes: consolidated materials (ROIs 4 and 5), talus or mass-wasted materials (ROIs 2 and 6), and smooth unconsolidated materials (ROIs 1 and 3). In Fig. 2 (lower panel) the location of ROIs is shown in an OSIRIS NAC image acquired of the same region with different illumination geometry to investigate the regional morphology on the surface.
The ROIs 4 and 5 are morphologically similar to each other at the investigated resolution, as expected because they represent rim materials for the same circular structure enclosing the area that contains ROIs 1 and 2. ROI 6 appears brighter in the low phase images, but there are no clear morphological differences. Topographically, ROI 2 materials reside on a lower slope than ROI 6 materials. Both are essentially dominated by boulders of various sizes and shapes. Finally, ROIs 1 and 3 are smooth at the investigated resolution and mainly comprise materials that dominate the interior of the Imhotep region. ROI 3 is predominantly a patch that appears to be brighter than its surrounding at higher phase angles ( Fig. 2 lower panel). ROI 1 is situated in a region where potential material loss occurred around perihelion (Groussin et al. 2015;El-Maarry et al. 2017), and it might therefore be composed of recently excavated or mobilized materials.

Surface phase curve
For each ROI we built the variation of surface reflectance with respect to a certain phase angle in three wavelengths, using the photometrically corrected images. Each pixel in the corrected images gives the dependence of the surface reflectance on phase angle (α), called the surface phase curve (Li et al. 2015). We here refer to the surface phase curve as the phase curve.
In order to associate the location of ROIs in each image sequence accurately with the corresponding phase angle, we converted the images and the phase angle distribution to a simple cylindrical projection with a resolution of 0.53 m pix −1 . The corresponding latitude and longitude coverage was calculated using the x,y, z coordinates from the global 3D shape model that was constructed based on the stereo-photogrammetry technique (Preusker et al. 2017).
The mean reflectance of each ROI was calculated from the map-projected RGB images, and the corresponding mean phase angle from the projected phase angle images was extracted. The resulting phase curves in three wavelengths for each ROI are plotted in Fig. 3. Each ROI contains 15 779 pixels. The uncertainty on the mean reflectance is represented by the standard deviation of the mean of data points within the ROI.

Linear-exponential modeling
To qualitatively classify and compare the phase curve properties in the small phase angle range, we made use of an empirical mathematical model, known as the linear-exponential model (Kaasalainen et al. 2001;Rosenbush et al. 2002;Muinonen et al. 2009). This model uses a four-parameter linear-exponential function to reproduce the phase curve. The function includes OE parameters, considering the phase curve as combination of an exponential peak and a linear part, and is given by where I/F s is the amplitude of the opposition peak and is defined as the reflectance increases relative to the background reflectance I/F b . B is the slope of linear part, and d is the angular width of OE. The angular HWHM of the OE is accordingly calculated as (Muinonen et al. 2002) The amplitude of the OE in the form of the enhancement factor, ζ is defined as (Muinonen et al. 2002;Rosenbush et al. 2002) A11, page 3 of 10 A&A 630, A11 (2019) Fig. 3. Phase curves of six ROIs that are extracted from the Ash-Khepry-Imhotep region of 67P at three wavelengths 480.7 nm (blue), 649.2 nm (green), and 743.7 nm (red). The phase curves span a small phase angle range (α < 10 • ). The dash lines represent the fitted linear-exponential model (Fig. 1). The dotted lines show a linear phase law.
We used a weighted nonlinear least-squares method in MATLAB to fit the function to the phase curves of each ROI for three wavelengths (see Fig. 3) to constrain the four parameters for each ROI and wavelength. The Levenberg-Marquardt algorithm (Levenberg 1944;Marquardt 1963;Moré 1977) was used for the fitting procedure. We further calculated the best-fit value A11, page 4 of 10 N. Masoumzadeh et al.: Phase-curve analysis of comet 67P/Churyumov-Gerasimenko at small phase angles of ζ according to Fig. 3. The model HWHM and ζ for each ROI and wavelength are listed in Table 1 and displayed in Fig. 4. The uncertainties for the four free parameters are the 1σ values returned by the fitting algorithm. The error for ζ is estimated based on the error propagation formula. Opposition effect. To clearly show the OE, which is defined as a departure of the phase curve from linearity toward zerophase angle, we show the linear part of the fitted linearexponential function by setting I/F s = 0. In Fig. 3 we plot the linear section that is I/F b + Bα. The departure from linearity is obvious for all derived phase curves and thus reveals the OE.
Spectral behavior of the opposition effect. The theory of CB predicts a variation in HWHM with wavelength, specifically, according to Mishchenko (1992), the HWHM (in degrees) changes with wavelength as where λ is the wavelength, Q sca is the scattering efficiency, f is the filling factor of the medium, and r o is the radius of particles. Because f and r o do not change with wavelength and the values of Q sca are very similar for the range of real (1.55-1.75) and imaginary parts (0.001-0.1) of the refractive index that cover most typical silicates and organics (Kolokolova et al. 2003), the HWHM mainly depends on the wavelength and increases as the wavelength increases. In order to determine whether the variation of the model OE parameters with wavelength is significant, we used a weighted linear fit to calculate the variation of the model HWHM and ζ versus wavelength for each ROI (see dashed lines in Fig. 4). The best-fit slopes for all ROIs are on the order of 10 −4 , while the resulting error on the slopes is estimated to be on the order of 10 −3 . We therefore infer from this analysis that there is no statistically significant evidence of a wavelength dependence on the OE characteristics.

Phase-ratio analysis
To estimate how much the phase curve at small phase angles depends on the structure and the roughness of the terrain, we studied the phase-ratio versus reflectance (Shkuratov et al. 2011) for each ROI. The phase ratio was constructed by calculating the ratio of the surface reflectance measured at two different phase angles: I/F(α 1 ) I/F(α 2 ) at α 1 ∼ 0 • and α 2 > α 1 . The phase ratio helps to suppress the albedo variations, leaving only changes related to phase angle. Using the phase-ratio technique, the influence of spatially unresolved roughness and microtopography on the phase-curve slope can be illustrated ).
An inverse correlation between phase ratio and reflectance values is expected for the SH, as discussed by Shkuratov et al. (2011). The anticorrelation appears because a surface with higher reflectance experiences enhanced multiple scattering that causes light to penetrate the shadows and thus decreases the shadowed area and the SH, which in turn causes a lower slope of the phase curve. Furthermore, a diagram that shows the anticorrelation between phase ratio and reflectance, as discussed by Shkuratov et al. (2012), can also distinguish between rougher and smoother terrains based on the deviation of the data from the orthogonal regression line. The larger deviation corresponds to the terrain with higher roughness.
We built phase-ratio images from the two map-projected images that were acquired at phase angles α 1 ∼ 0 • and α 2 ∼ 5 • together with the corresponding reflectance map at larger phase angle α 2 for each ROI, as illustrated in Fig. 5. Darker areas are linked to the shallower phase curve and brighter areas correlate with higher values of the ratio, which are related to the steeper phase curve (Shkuratov et al. 2011;Kaydash et al. 2012).
In Fig. 6 the phase-ratios versus reflectance of the surface at angle α 2 for each ROI are illustrated by a 2D histogram. We note that the reflectance at α 2 ∼ 5 • that is lower than 0.04 was excluded in our analysis to avoid shadowed regions that could lead to an unphysically high phase-ratio value. The 2D histogram uses the bivariate normal probability distribution (Wilks 2011) to group data into the 2D bins, and each bin is colored based on the frequency of occurrence. The binning method is based on the Scott rule (Scott 2010), which uses a bin size of 3.5 × σ(x) x 1/4 , 3.5 × σ(y) y 1/4 , where σ is the standard deviation. The bivariate normal distribution is able to describe the overall shape of data with properties that depend on the mean vector and the covariance matrix of the two data sets. The mean vector corresponds to the coordinate of the most frequent occurrence in the two groups of data. The spread and orientation of the scatter data can be characterized by the covariance matrix. In order to visualize the distribution of phase-ratio reflectance histograms, we computed a 95% confidence ellipse derived from the covariance matrix and centered on the mean vector of the phase-ratio data set and the reflectance data set A11, page 6 of 10 N. Masoumzadeh et al.: Phase-curve analysis of comet 67P/Churyumov-Gerasimenko at small phase angles for each ROI. The confidence ellipse (red solid line in Fig. 6) has two properties of the width and the orientation that provide information about the pattern and the correlation between the phase ratio and the reflectance. The oriented confidence ellipses with a slope value of −88 • for all ROIs indicate an inverse correlation between the phase ratio and reflectance. A11, page 7 of 10 A&A 630, A11 (2019) Fig. 6. 2D histogram of phase ratio vs. reflectance at angle α 2 for each ROI at λ = 649.2 nm. The phase ratios were collected from images acquired at α 1 ∼ 0 • and α 2 ∼ 5 • . The corresponding reflectance was measured for the higher phase angle, α 2 . The color scale represents the number of data points within the bin, from high (yellow) to low (dark blue). The solid red line shows the 95% confidence ellipse.
This comes from the fact that the covariance matrix is not diagonal.
Visual inspection of the confidence ellipses in Fig. 6 suggests that the phase-ratio data points of ROIs 1 and 3 with a narrow reflectance range fill the confidence ellipse and shows less spread outside the ellipse. ROIs 2 and 6 display a wide reflectance range and larger scatter. The data points in ROI 6 include higher phase ratios than ROI 2, suggesting a rougher terrain for ROI 6. This is consistent with our description in Sect. 2.2, where we indicate that the surfaces of ROIs 1 and 3 are rather smooth, but ROIs 2 and 6 are talus regions that represent a rougher surface.
Although ROIs 4 and 5 are both classified as consolidated regions, differences between the spatial distributions in their phase-ratio reflectance histograms are observed. As illustrated in Fig. 6 for ROI 5, the phase ratios corresponding to the lower reflectance values are more scattered and deviate from the confidence ellipse, while the phase ratios of higher reflectance are clustered inside the confidence ellipse. We interpret this behavior in phase-ratio versus reflectance histogram of ROI 5 as caused by differences in the surface structure of ROI 5, for instance, subresolution roughness or differences in grain size distribution, which result in the differences in scattering behavior.
A special spectral behavior for ROI 5 was also noted by Feller et al. (2019). ROI 5 corresponds to the Cuesta feature that was analyzed by these authors. Using the wavelength range 535-743 nm, they found a low spectral slope and higher reflectance for the Cuesta feature than for its surrounding. They suggested that this spectral behavior is evidence of distinct compositional properties.

Summary and conclusion
We analyzed the OE behavior for the nucleus of comet 67P. We extracted the phase curves in the small phase angle domain (0 • < α < 10 • ) for six regions of interest. The phase curves were built in three wavelengths from OSIRIS NAC images that cover the Ash-Khepry-Imhotep region.
We found a departure from a linear phase law that accounts for the OE. We qualitatively studied the OE characteristics with respect to the wavelength using a linear-exponential function. No strong wavelength dependency is observed in all best-fit OE parameters for all ROIs.
We applied a phase-ratio analysis to the ROIs to study the structural properties of the surface. An inverse correlation is observed between the phase ratio and reflectance for ROIs. This means that a lower phase-curve slope is typical for brighter materials because the shadowed area is reduced as a result of multiple light scattering. This implies that the main effect that defines any opposition effect of comet 67P is the SH.
This anticorrelation behavior is also observed in lunar images . No sign of any correlation between the phase ratio and reflectance for mercurian surfaces was found by Blewett et al. (2014). Although factors other than composition, such as surface structure, might cause the effect of multiple scattering to weaken, the authors speculated that the low reflectance of materials on the mercurian regolith is responsible in this case.
We found a good agreement between the three morphological classes of the defined ROIs on the Ash-Khepry-Imhotep region when we plotted a phase-ratio reflectance histogram. We argue that the scatter outside of the confidence ellipse in the phase-ratio reflectance histogram of ROI 5, a consolidated region, may be connected to the structural and compositional properties of the region.
The full understanding of the opposition effect requires a robust light-scattering model and combination of the information from different instruments, including ground-based observations (Snodgrass et al. 2011;Kokotanekova et al. 2017). Several computational techniques are available (Mackowski & Mishchenko 1996, 2011Muinonen et al. 2012Muinonen et al. , 2018) that suggest different physical and mathematical complexities that need to be overcome to approach this problem. The data we presented in the paper and the qualitative analysis we performed not only allowed us to suggest some conclusions about physical and compositional characteristics of the nucleus of comet 67P, but also to represent a unique set of small-phase angle data that can be used by the modelers to test and improve the validity of their computational approaches and to compare different techniques to solve the light-scattering inverse problem.