The opposition effect of 67P/Churyumov–Gerasimenko on post-perihelion Rosetta images

Abstract

High-resolution OSIRIS/Rosetta images of 67P/Churyumov–Gerasimenko acquired on the night run of 2016 April 9–10 show, at large scale, an opposition effect (OE) spot sweeping across Imhotep as the phase angle ranges from 0° to 17°. In this work, we fitted the phase curve of the whole surface imaged as well as three particular features using both the linear–exponential and Hapke models. These features encompass different types of spectral behaviour: a circular mesa, one venous structure and an assemblage of bright spots, going from red to blue colours. Both the Hapke and linear–exponential parameters indicate a stepwise sharpening of the OE from bright spots to circular mesa. Yet a very broad nonlinear phase curve is verified and no sign of sharp OE associated with a coherent-backscattering mechanism is observed. We estimate that the 67P surface is dominated by opaque, desiccated and larger-than-wavelength irregular grains. Veins and bright spots display photometric properties consistent with surfaces becoming slightly brighter as they are enriched by high-albedo ice grains. We also report the estimation of normal albedo for all cometary regions observed throughout the image sequence. Comparison to pre-perihelion results indicates that far better insolation of northern brighter regions, i.e. Hapi, Hathor and Seth, is sufficient to explain mismatches on the photometric parameters. However, metre-scale photometric analysis of the Imhotep–Ash boundary area advocates for mild darkening (<7 per cent) of the surface at local scale.

Rosetta: 67P/Churyumov–Gerasimenko, photometry

1 INTRODUCTION

The Rosetta mission came to an end on 2016 September 30 after a fruitful study of more than two years of the comet 67P/Churyumov–Gerasimenko. The rendezvous took place in 2014 August, when Rosetta concluded the approaching phase and progressively adjusted its orbit around the comet. The surface of the nucleus was surveyed by the Optical, Spectroscopic and Infrared Remote Imaging System (OSIRIS), the ‘eyes of Rosetta’, composed of two cameras (NAC, narrow angle camera and WAC, wide angle camera) with 21 broad- and narrow-band filters ranging from 269 to 989 nm (NAC) and from 246 to 629 nm (WAC). The OSIRIS imaging system is described in length by Keller et al. (2007).

Rosetta had three opportunities to capture images of the nucleus at very small phase angles: on 2014 July 29, though images taken along the approach (Fornasier et al. 2015); on 2015 February 14, during the pre-perihelion fly-by (Feller et al. 2016; Masoumzadeh et al. 2017); and recently, on 2016 April 10, during the distancing stage to a second fly-by. Observations at small phase angle had been generally avoided for security measures. Opposition configuration inevitably leads to 180 degrees of phase angle to achieve an orbital closure that could temporarily shut the line of communication with Rosetta’s antenna and leave the spacecraft in the dark for several minutes. Fortunately, Rosetta successfully performed all recovery manoeuvres and has delivered important phase-curve measurements at very small phase angles that hint at the uppermost layer properties.

The so-called photometric phase curve is a function of the reflectance with respect to the phase angle between the surface-observer vector and the surface-source vector. The main morphological aspect of a phase curve of compact random discrete media (i.e. dust cover, regolith, sand and so forth) is the opposition effect (OE), a nonlinear increase in brightness when the phase angle approaches zero degrees. It is connected to two convolved optical mechanisms: the shadow-hiding effect (SHOE, e.g. Hapke 1981; Shkuratov 1983; Lumme, Peltoniemi & Irvine 1990; Stankevich, Shkuratov & Muinonen 1999) and coherent-backscattering enhancement (CBOE, e.g. Albada & Lagendijk 1985; Akkermans, Wolf & Maynard 1986; Muinonen 1994; Mishchenko et al. 2009). In the former mechanism, shadows amid larger-than-wavelength opaque particles are progressively hidden to the observer as the geometry approaches opposition. In the latter, the multiple-order scattered electromagnetic waves of a cluster of smaller-than-wavelength scatterers constructively interfere at very small phase angles. These mechanisms are intrinsically connected to optical indexes, packing factor, size distribution, particle shape, inclusions, chemical composition and opacity/transparency in a complex interplay that is still the subject of ongoing research (e.g. Shkuratov et al. 2002; Nelson et al. 2002; Kaasalainen 2003; Muinonen et al. 2012; Déau et al. 2013b). A secondary major aspect is the macroscopic shadowing, a mechanism often significant when surfaces are observed at phase angles larger than ∼60°. Boulders, micro-craters, protuberances and micro-irregularities are generally evoked to explain the hindering of brightness due to castings of long shadows (Hapke 1984; Buratti & Veverka 1985; Goguen et al. 2010; Shkuratov et al. 2012). The phase curve also depends on the incidence and emergence angles, as they are closely related to particle scattering behaviour and, macroscopically, to the extension of observed shadows.

The comet 67P/Churyumov–Gerasimenko presents a dark nucleus with geometric albedo of 6.5 ± 0.2 per cent at 649 nm (Fornasier et al. 2015). The linear slope of the 67P photometric phase curve at 3.3 AU (Fornasier et al. 2015) is very similar to those found for Tempel 1 (Li et al. 2007a, 2013), Hartley 2 (Li et al. 2013) and Borelly (Li et al. 2007b) and close to those from low-albedo asteroids and Martian satellites (Masoumzadeh et al. 2017). High steepness on the phase slope and absence of sharp OE is also observed among some dark objects of the Solar system (Shevchenko et al. 2008, 2012) and is generally interpreted as prevailing single-scattering effects. Additionally, previous Hapke modelling of OSIRIS images show a very broad SHOE-H,¹ indicating that 67P dust grains are possibly sufficiently opaque to disregard the strong contribution from inter-grain multiple scattering and the coherent-backscattering mechanism (Fornasier et al. 2015; Ciarniello et al. 2015; Feller et al. 2016). Results from RADAR techniques also point to a dry and highly porous dust cover composed of irregular opaque micrometric particles (Kamoun et al. 2014; Lethuillier et al. 2016). Recently, individual cometary grains of tens to several hundred micrometres collected by the COSISCOPE (spatial resolution of 14 μm) instrument showed unambiguous intra-grain shadow casting when illuminated by an LED at a grazing incidence angle with opaque granular structures at a scale of about 30–40 times larger than the wavelength (Langevin et al. 2016; Hilchenbach et al. 2016).

In this work, we report the analysis of high-resolution images exhibiting the OE phenomenon and some peculiar features found therein. Never before has the full opposition spot been pictured on a global scale among cometary nuclei. ‘Phase zero’ images are glimpses of the true albedo (Shkuratov et al. 2011) of a cometary surface and are directly correlated to composition. Therefore, our goal is to understand the radiative scattering properties of the ice-rich features and of the whole comet at the extent of the regions observed therein. By characterizing and retrieving meaningful photometric parameters, we can help to constrain the amount of exposed ice on the comet and the properties of the dehydrated and organic-rich material that coats the nucleus.

2 OBSERVATIONS AND AUXILIARY DATA

A sequence of images exhibiting the OE spot was acquired between 2016 April 9 and 10 with a WAC F18 filter of the OSIRIS camera (central wavelength of 612.6 nm and bandwidth of 9.8 nm). The full set comprises 10 images spanning 3 h and 8 min, whereas 7 images with a duration of about 1 h show the small-phase-angle event. Fig. 1 shows the images that exhibit OE along with the projection of the phase angle over the same images. The spacecraft (hereafter referred as S/C) distance and resolution ranged from 27.5 km (2.8 m px⁻¹) to 30.23 km (3 m px⁻¹). The images have their line of sight starting under Apis–Atum, sweeping over Imhotep–Khepry–Ash until it reaches the Khepry–Aker–Bastet boundaries. The opposition surge spot was visible under Imhotep–Khepry, whereas the last two images spotted the Bastet–Hathor interface. Fig. 2 shows the boundaries of the regions in the field-of-view at mid-opposition surge transit. Further morphological descriptions of the regions on the 67P nucleus are given in Thomas et al. (2015), El-Maarry et al. (2015) and El-Maarry et al. (2016).

Figure 1.

Morphological regions and their boundaries of various colours. The figure represents all regions observed in the field-of-view of 2016 April 9 UT 23:59:32, the time date of the central image to the OE event.

Figure 2.

WAC F18 image sequences from 2016 April 9 UT 23:29:32 to 2016 April 10 UT 00:32:32 exhibiting the OE phenomenon with their corresponding phase-angle graphical representations. Top seven grey images: calibrated photometrically uncorrected RADF images. Bottom seven heat-scale images: phase-angle gradient in degrees. The OE spot starts at Apis and sweeps over Imhotep until reaching Bastet and finally leaving the surface.

To better constrain the morphology of the phase curve, we added to our sample 8 images acquired through the WAC F18 filter exactly 2 months before, on 2016 February 10. Similar to the 2016 April 10 sequence, Imhotep and neighbouring regions are on the field-of-view at 7 images. Their phase angle covered 62.0° to 67.9° over 3 h and 35 min. Their S/C distance and resolution varied from 48.5 km (4.9 m px⁻¹) to 48.1 km (4.85 m px⁻¹) as Rosetta flew from right above Imhotep to Anhur–Khepry–Sobek, at the Southern hemisphere. Table 1 contains information on the sub-S/C coordinates, phase angle, phase-angle range, spatial resolution, solar distance and S/C distance of the images of both sequences.

Table 1.

Characteristics of the WAC F18 images from the 2016 February 10 and April 10 sequences.

Time	αsub-S/C (°)	Δα(°)	D (km)	r (AU)	Res. (m px−1)	Sub-S/C latitude	Sub-S/C longitude
2016-04-09T21:47:33.392	15.12	10.3–17.5	27.75	2.766	2.8	−8°31΄12″	28°33΄36″
2016-04-09T23:29:32.760	2.54	0.0–6.3	27.43	2.767	2.78	−8°23΄24″	−6°43΄48″
2016-04-09T23:38:32.750	1.56	0.0–5.4	27.49	”	2.78	−8°24΄36″	−9°55΄48″
2016-04-09T23:50:32.738	0.65	0.0–4.7	27.60	”	2.79	−8°26΄24″	−14°12΄00″
2016-04-09T23:59:32.754	1.18	0.0–5.8	27.69	”	2.80	−8°27΄00″	−17°24΄36″
2016-04-10T00:11:32.773	2.44	0.0–7.3	27.83	”	2.81	−8°28΄12″	−21°41΄24″
2016-04-10T00:20:32.738	3.43	0.0–8.3	27.95	”	2.82	−8°28΄48″	−24°54΄36″
2016-04-10T00:32:32.993	4.75	0.7–9.5	28.12	”	2.84	−8°29΄24″	−29°12΄36″
2016-04-10T01:28:40.738	10.82	6.1–14.8	29.16	”	2.94	−8°28΄48″	−49°32΄24″
2016-04-10T01:55:49.766	13.63	9.1–17.3	29.76	”	3.01	−8°26΄24″	−59°30΄36″
2016-02-10T13:53:51.616	65.00	63.1–67.8	47.88	2.331	4.84	27°0΄0″	−93°54΄0″
2016-02-10T14:01:54.768	64.99	63.1–67.8	47.86	2.332	4.83	26°55΄12″	−90°6΄36″
2016-02-10T15:20:18.543	64.85	63.9–67.8	47.54	”	4.80	26°8΄24″	53°10΄48″
2016-02-10T15:28:17.860	64.84	63.9–67.7	47.51	”	4.80	26°3΄36″	49°25΄12″
2016-02-10T16:20:18.604	64.76	62.4–67.2	47.22	”	4.77	25°31΄48″	24°55΄12″
2016-02-10T16:28:17.796	64.74	62.5–67.1	47.19	”	4.77	25°27΄36″	21°9΄36″
2016-02-10T17:20:18.568	64.66	62.6–66.7	46.99	”	4.75	24°55΄48″	−3°20΄24″
2016-02-10T17:28:19.065	64.65	62.7–66.7	46.97	2.333	4.74	24°51΄0″	−7°6΄36″

Relative errors for the absolute coefficient of WAC F18 and other filters were estimated by Tubiana et al. (2015). The absolute coefficient is a factor multiplied by DNS values to obtain the intensity and main source of uncertainties in the radiometric data. This way, WAC F18 is affected by an error no higher than 0.6 per cent. Photometric correction increases the uncertainties by accumulating errors associated with S/C pointing and mostly to shape-model approximations. This subject is further explored in Appendix C. All images were then corrected for camera optical distortions, radiometrically calibrated 32-bit float data and converted to radiance factor (r_F or RADF, dimensionless) units (Tubiana et al. 2015).

High-resolution shape models of the target, S/C orientation, target and S/C trajectories, instrument set-ups and rotational axis of the target are all fundamental to retrieve the incidence (i), emergence (e), azimuth (φ) and phase angles (α) for a detailed photometric analysis. We used the shape model SHAP 8 version 1.8 provided by Jorda et al. (2016) and developed through the stereophotoclinometry technique (SPC). SHAP 8 1.8 is presented as a raster digital elevation model of 3.14 million facets with a lateral resolution of 6 m per facet. The spherical coordinate system is referenced to the Cheops boulder on Imhotep. We avoided cliffs by removing facets at the shadows’ interfaces. The trajectories and instrumental information for Rosetta and 67P/Churyumov–Gerasimenko used to compute the state of the comet on each image are available through NAIF SPICE kernels² (Acton 1996) maintained by European Space Agency (ESA). The OASIS image simulator (Jorda et al. 2010) was used to reproduce the OSIRIS images and retrieve the photometric angles per shape-model facet. Among its outputs, OASIS provides rendered FITS images on the same observational conditions as the original OSIRIS images. We use these rendered images as reference to inconsistencies regarding any translation compared to the original images. We registered a consistent offset and rotation of X = −68.8 ± 2.3, Y = 20.8 ± 14.2 and R = −0.01° ± 2°. Thereunto, we translated the OSIRIS images to the OASIS images using the discrete Fourier transform³ (Reddy & Chatterji 1996). The co-registered images then gave us the best overlap between the OSIRIS and OASIS pixels for each image and therefore the best pixel–facet link (Hasselmann et al. 2016).

At very small phase angles the Sun must be counted as an extended source due to the non-parallelism of the incident rays. At this level, what was once approximately considered a point source now must be treated as diffusive. The angular size of the Sun imposes a lower limit to the phase angles that we can sample (Shkuratov 1991). Such a limit is trivially calculated by trigonometry and is given by |$\arcsin (\mathrm{R}_{\odot }/r)$|⁠, where R_⊙ is the radius of the source (R_⊙ = 6.957 × 10⁸ m) and r is the distance from object to source. Consequently, we removed all measurements at α ≲ 0.095°.

3 REGIONS OF INTEREST

As the OE spot sweeps across the surface, we observe several local albedo variations and some peculiar features. To study the phase curve of specific areas in the nucleus, we selected three regions of interest (ROI) based on their morphology or contrast in albedo that are observed throughout both sequences. In Fig. 3, where the image 2016 April 9 UT 23:59:32 is shown as an example, we identify each ROI and many other albedo variations, mainly dark rough terrains and smooth dusty seas like the one in central Imhotep. In Fig. 4, we present in false RGB colour a partial view of the Imhotep region and its corresponding spectral slope map based on images obtained with the OSIRIS NAC during the same 2016 April 10 observations. We can discriminate that blue spectral features are also the brighter ones. Imhotep has mostly red-to-average spectral slope, while two of the selected ROI are connected to the blue spectral units (Fornasier et al. 2015; Oklay et al. 2016b). The ROI are described in detail in the subsections below. They are defined and labelled according to previous works on morphological features in the Imhotep region: El-Maarry et al. (2015), El-Maarry et al. (2016) and Auger et al. (2015).

Figure 3.

(a) Contrast-stretched image from 2016 April 9 UT 23:59:32. The large central bright spot on the surface corresponds to the enhancement due to small phase angles. ROI 1, 2 and 3 are marked by colour boxes. (b) The estimated normal albedo for the same image, as a percentage. Topographic–photometric correction using the Hapke IMSA model (see Appendix B). Measured RADF are normalized by the radiance factor calculated from the model. All radiative dependences related to i, e and α are removed to an uncertainty of about 4 per cent.

Figure 4.

Composition of NAC 2016 February 10 UT 15:20 and 15:28 colour sequences revealing the Imhotep region and the regions of interest of this study. The sequence has spatial resolution of 0.9 m px⁻¹ and central α = 64.9°. (a) RGB image produced by the STIFF software using the F24 (480 nm), F83 (535.7 nm) and F41 (882.1 nm) filters. Bluer colours represent lower spectral slope. (b) Spectral slope maps of the same colour sequence. The images are topographically corrected using the Lommel–Seeliger law and co-registered using projective homography and dense optical flow (Farnebäck 2003). The spectral slopes are calculated using the equation |$S[{\rm {per \,cent}}/100\ \mathrm{nm}]=\frac{R_{882}-R_{535}}{R_{535}}\cdot \frac{10^{4}}{(882-535.7)}$| (Fornasier et al. 2015).

3.1 ROI 1: Circular Mesa

The Circular Mesa is the largest circular feature on Imhotep (latitude = −35°, longitude = 165°55^΄48^″), extending for about 647 metres of diameter (El-Maarry et al. 2015; Basin F in Auger et al. 2015). The structure is partially covered by the very same granular material that comprises most of Imhotep dusty deposit (see Fig. 4a and Fig. 3). On the opposite higher wall, the Mesa exhibits an irregular fractured pattern and some scattered debris, demonstrating that the feature has undergone some landslides in the past. The Mesa also harbours a small area of higher albedo that extends for 47.6 metres in diameter (6 per cent of total surface) with an aspect very similar to that of the blue veins (ROI 2; see below). Auger et al. (2015) proposed that the large circular features are either impact craters or the rising up of gas/void chambers from the interior of the nucleus that emerge from the thinning of the upper layers. The spectral slope map (Fig. 4b) shows that the Circular Mesa has a spectral slope of ∼18 per cent/100 nm at α ≈ 65° in line with Imhotep as a whole. The small bright patch, on the other hand, has a lower spectral slope (S ≈ 15 per cent/100 nm).

3.2 ROI 2: blue veins

This venous bright area of mirror-like appearance was spotted throughout the OE images (see Figs 1 and 3). Located in the vicinity of the Circular Mesa, the blue veins resemble exposed veins of an ice-rich layer (−15°30^΄, 156°30^΄). The area agglomerates several small roundish features that may be related to the origin of this patch (see Fig. 4a). Auger et al. (2015) proposed a scenario where such small roundish features were ancient degassing conduits that were further exposed through the progressive erosion of the surrounding surface. Richer in ices, this eroded layer would be expected to be bluer in colour (Oklay et al. 2016b; Barucci et al. 2016). Indeed, this is actually observed through spectrophotometry (S ≈ 14.5 per cent/100 nm); see Fig. 4b). Recently, Oklay et al. (2016b) and Oklay et al. (2016a) have worked on the spectrophotometry of the same area, which was revealed to be persistently bluer than average since 2014 September. Furthermore, Knollenberg et al. (2016) have demonstrated that the blue veins were also the source region of an outburst detected on the Ides of March 2015. Therefore, ROI 2 can be considered an active area up to about one year before our observations.

With respect to its apparent albedo, the blue veins are brighter than their vicinities for α ≲ 10° but fade into the same average photometric appearance of Imhotep at α ≈ 65°.

3.3 ROI 3: bright spots and patches

Bright spots or patches, found throughout several regions of 67P, have been largely studied and documented (Pommerol et al. 2015b; Barucci et al. 2016; Oklay et al. 2016b; Deshapriya et al. 2016; Fornasier et al. 2017). Such features of high radiance-factor contrast can be spotted from a few pixels wide to a full identifiable structure, depending on the spatial resolution and observational conditions of the images (Barucci et al. 2016; Filacchione et al. 2016). A bright feature is generally located bordering scarps, cliffs or part of some crumbling material. Their apparent albedo can be to 10 times higher than the neighbouring pixels at larger phase angles. These features are thought to be either remnants of exposed H₂O ice-rich layers, partially exposed underground ice-rich layers or ephemeral H₂O frost replenishment (Fornasier et al. 2016), depending on their surviving time, morphology and their behaviour according to insolation and environmental shadow casting. Barucci et al. (2016) studied some of their available VIRTIS spectra (Coradini et al. 2007) and found that water-ice abundance is about a few per cent. Fornasier et al. (2016) state that water ice may be up to 30 per cent for one specific spot. With respect to CO₂ signatures, they are mostly difficult to identify on the surface due to its lower sublimation point compared to H₂O ices. Nonetheless, Filacchione et al. (2016) recently studied one ice-rich area in the Anhur region and found 0.1 per cent of CO₂, the first detection of CO₂ ice on a cometary surface.

We then selected 78 unique spots available in our observations to compose an average phase curve. Gathering and analysing all bright spots together is necessary for having a satisfactory (i, e, α) coverage. Fig. 5(b) shows latitude and longitude distributions of the bright spots. The largest bright patch was located on the Ash region (−47°26^΄24^″, 91°3^΄) and is the single identifiable bright feature throughout all the 2016 April 10 sequence. This patch is about 147 metres in diameter and is a cluster of 8 bright spots that were individually selected to make up our sample. Most of the spots and patches were observed during the full extent of one hour in the OE sequence and seemed stable during this time range. The few disappearances are well correlated to slight changes of the field-of-view and their hiding behind cliffs. Two features bearing r_F > 0.10 are observed once and may be related to the short timescale of frost deposition.

Figure 5.

ROIs 2 and 3: Blue veins and bright spots. (a) Blue veins as depicted in the 2016 April 10 sequence. The degassing conduit features and average albedo terrain are filtered (see subsection 5.2). For the images taken at 2016 April 10 UT 01:28:40 and UT 01:55:49, the observational conditions were too oblique (e > 70°). Grey-scale in RADF. (a) Longitude and latitude of gathered bright spots. Colour bar of the figure on the left represents the maximum radiance factor of the spot, whereas on the right we have the phase angle on which the spot was recorded.

4 METHODOLOGY

The methods that we used to give an interpretation to the phase curves have been extensively reported in the literature. Our approach is to undertake two different analyses and give two different sets of parameters for describing our data. First, we fit the phase curve using a simple linear–exponential formula of four free parameters (Kaasalainen, Muinonen & Piironen 2001; Muinonen et al. 2002; Kaasalainen et al. 2003; Rosenbush et al. 2005) after correcting the (i, e) angles by the Lommel–Seeliger law: the amplitude of the opposition surge A, the width of the opposition surge d, the phase curve albedo without opposition surge b and the angular coefficient of the linear part k. We apply the linear–exponential formalism to retrieve parameters based only on the morphology of the phase curve. In Appendix A, we outline the formula and the minimization technique for estimating parameters. Secondly, we employed the most recent version of the Hapke isotropic multiple-scattering approximation (IMSA) model (Hapke 2012) to estimate the collective characteristics of the optically active layer of the nucleus. Hapke (2012) explicitly incorporates the role of superficial porosity (1 − ϕ) into the model, improving the estimation of the single-scattering albedo w, asymmetric factor g_sca, SHOE-H width h_SH and SHOE-H amplitude B_SH. The Hapke IMSA model is outlined in Appendix B. The way in which uncertainties are calculated for both approaches is described in Appendix A. We have also conducted a detailed analysis of errors associated with topographic–photometric correction by the Lommel–Seeliger law in Appendix C.

5 RESULTS

5.1 Phase-curve analysis

The complete set of images constitutes a full phase-curve stretching to 0.095° < α < 67.9° with a gap at 17.1° < α < 62.0°. A total of 997 701 facets of SPC SHAP 8 were observed. To reconstruct the phase curve for each ROI, we manually selected all their pixels and corresponding facets. For the Circular Mesa, we conserved all facets, whilst for the blue veins and the bright spots we filtered every facet within a 1σ envelope surrounding the global solution (σ_RADF = ±0.0025). This was arbitrarily done because our scope of work is merely to investigate the spectrally neutral and brighter components, as illustrated in Fig. 5(a). The bright spots are too small to be reproduced in a similar figure. In particular, we decided not to apply the Lommel–Seeliger disc correction for the bright spots since most of the features are located close to cliffs or sharp boundaries and extend for just a few pixels. Therefore, accumulative errors associated with topographical photometric correction could unsatisfactorily reach higher than 10 per cent.

Subsequently, we binned the (i_n, e_n, α_n) table of n facets into a 20 × 20 × 50 cell grid for each image, corresponding to 3.7° × 3.7° × (0.09° − 0.35°) steps. The radiance factor is averaged for each cell, mitigating any influence of variegation and poor pixels/facets (Hasselmann et al. 2016). Once more, we suppressed this step for bright spots due to discontinuities in the gathered facets. The results from modelling through the linear–exponential equation and IMSA are presented in Table 2.

Table 2.

Linear–exponential and porosity-dependent Hapke IMSA parameters. The symbols are described in Appendices A and B. |$\bar{\chi }^{2}$| is the average chi-squared and N represents the total number of facets.

Parameters	σ	All data	Mesa	Blue veins	Bright spots
A	±0.0008	0.0377	0.0381	0.0305	0.0278
b	±0.001	0.0240	0.0233	0.036	0.0423
d	±0.004	0.172	0.184	0.152	0.108
k	±0.001	0.017	0.015	0.025	0.026
ρν, l–e ( per cent)	±0.05	6.17	6.14	6.63	7.01
Ia	±0.09	2.57	2.65	1.86	1.67
HWHM (°)	±0.23	6.79	7.27*	6.00	4.27
|$\bar{\chi }_\textrm{l--e}^{2}(\times 10^{-4})^{b}$|		6.4	2.6	6.5	3.1
w	±0.001	0.027	0.033	0.035	0.047
B0	±0.3	2.42	2.41	2.63	2.38
hs	±0.005	0.081	0.072	0.079	0.06
gsca	±0.02	−0.424	−0.38	−0.368	−0.335
|$\bar{\theta }$| (°)	±3°	26	21	33	[15]
ρν, Hapke ( per cent)	±0.05	6.14	6.23	6.66	7.27
K	±0.005	1.245	1.234	1.238	1.198
1 − ϕc ( per cent)	±2 per cent	82	84	83	86
HWHMSH (°)	±0.3°	9.28	8.82*	8.94	6.88
|$\bar{\chi }_\textrm{l--e}^{2}(\times 10^{-4})$|		8.7	3	3.4	4.1
N		997 701	27 209	85 174	28 879

|$^{a}\ I=\frac{(A+b)}{b}$| is the amplitude of the exponential term.

|$^{b}\frac{1}{N}\sum\limits_{i=0}^{N}(r_{F,i}-r_{\mathrm{model},i})^{2}/r_{\mathrm{model},i}$|⁠.

^c ϕ represents the filling factor of the optically active layer; thus 1 − ϕ is then called superficial porosity and is not related to the porosity of the whole nucleus.

In Fig. 6, we present the phase curves for all data and each ROI alongside their respective photometric modelling. The trend on the χ² residuals as a function of phase angle is no larger than 5 per cent for every analysis. Inspection of the data at α < 5° shows a monotonic increase in the phase curves when disregarding albedo variation due to different terrains. No consistent sharp OE at α < 3° associated with the coherent-backscattering mechanism is observed.

Figure 6.

Phase-curve modelling using the Hapke IMSA model and linear–exponential fitting with Lommel–Seeliger disc correction. (Column a) Phase curves. Measured RADF are displayed either in cells or facet points. RADF computed through the Hapke IMSA model are shown by red dots, whilst the best linear–exponential solution is represented by a traced line. χ² residuals for both methods are shifted downward to −0.015 (IMSA) and −0.025 (L–E). Specifically, for the Circular Mesa, the mismatch between the photometric disc correction obtained from the Hapke IMSA and Lommel–Seeliger leads to differences in the estimated HWHM. There, green dots represent A_eq from the Lommel–Seeliger correction. (Column b) Phase curves at α < 5°. Colour labels similar to column (a). Linear–exponential fitting leads to a slight sub-estimation of normal albedo for ROI 2 and ROI 3.

Particularly with regard to the Circular Mesa, we observe a mismatch in the photometric correction for data taken at oblique emergence angles (e > 70°) using the Hapke IMSA model and Lommel–Seeliger law. This leads to consistent ±0.5° shifts in the half-width at half-maximum (HWHM) (marked with * on Table 2) though it remains close to the HWHM_all-data.

The photometric phase curves of the bright spots display far larger dispersion in RADF. As these features are mostly a few pixels wide under the spatial resolution of the WAC, the state of evolution and desiccation for each spot is hardly accessible. We consider that putting together mixed spots of diverse dust/ice ratios is responsible for such dispersion.

5.2 Morphological regions and other features

Photometrically corrected images enable us to identify albedo variations and correlate them with the morphological regions. Hence we used the Hapke IMSA model (Appendix B) to correct all pixels under i < 75° and not hidden in shadows. We have made this analysis only over the 2016 April 10 sequence. Fig. 3 corresponds to the best example of photometric correction, since the field-of-view does not vary greatly throughout the opposition images.

The dark rough terrains are mostly located at Imhotep and Bes. Wide rough terrains are generally 4–5 per cent darker and harbour few pixel-wide boulders closely packed together. We see two large portions of such terrains; one is close to the blue veins (3°47^΄, 172°58^΄) and the other is situated on the Imhotep–Bes boundary (20°38^΄, 124°16^΄) at the start of the Aten terraces (El-Maarry et al. 2016). We also identify a very large and dark boulder (−19°36^΄, 98°36^΄) about 90 metres wide on the Imhotep–Ash boundary and containing one bright spot (+30 per cent) on its border. Other dark structures are scattered in Ash and Bes, closer to the image limbs, and are likely shadowed terrains. All these features are visible on Fig. 3 (use the map on Fig. 2 as reference).

On later images (7° < α < 14°), Hapi, Hathor, Babi, Ma’at and Bastet become partially visible. We also have a clear view of Aten’s smooth terraces. The smooth terrain at Hapi is up to +8.8 per cent brighter than average, the same as the highest values found on the terraces. Ma’at is also brighter than average, but to a lesser extent (+4 per cent). A list of estimated median r_F/r_Hapke for each visible region, together with their 1st and 3rd quartiles, is shown in Table 3. We selected a box contained within the morphological regions and retrieved all active facets within. We remind the reader that not all regions are observed at very small phase angles (all except Apis, Imhotep, Khepry, Aten and Bastet), so they may have their r_F/r_Hapke underestimated by an unknown few per cent due to local shadows not accounted for by the shape-model resolution. To obtain the normal albedos, the ratios have to be multiplied by the all-data value ρ_ν,Hapke of Table 2.

Table 3.

Observed morphological regions and their corresponding albedo ratios. (lat, long)_central and (lat, long)_box represent the central position and width of the box from which N_facets was gathered and the median r_F/r_Hapke was calculated. w/w₅₅₀ is the relative single-scattering albedo from Filacchione et al. (2016) with respect to w₅₅₀ = 0.052 estimated by Ciarniello et al. (2015). Terrain classification is reproduced from Filacchione et al. (2016): SC = strongly consolidated; NCS = non-consolidated smooth; NCD = non-consolidated dust; DCB = dust-covered brittle; WCB = weakly consolidated brittle; D = depression.

Region	Terrain type	latcentral	longcentral	(lat, long)box	Nfacets	Nobs	rF/rHapke	w/w550
Anubis	NCS	−5°36΄36″	−120°57΄36″	38°18΄	21,320	1	|$0.951_{0.817}^{1.066}$|	0.83 ± 0.33
Anuket	SC	6°4΄12″	−61°49΄12″	50°58΄	117,168	1	|$0.976_{0.919}^{1.027}$|	0.87 ± 0.27
Apis	SC	4°1΄12″	−161°45΄0″	24°49΄	28,696	5	|$0.981_{0.959}^{1.001}$|	0.96 ± 0.31
Ash	NCD	29°6΄36″	154°52΄12″	60°39΄	139,528	9	|$1.006_{0.974}^{1.026}$|	0.94 ± 0.40
Aten	D	27°55΄48″	107°25΄48″	58°40΄	86,464	9	|$1.017_{0.999}^{1.037}$|	0.98 ± 0.31
Atum	SC	−6°21΄36″	−143°41΄60″	44°55΄	54,583	1	|$0.944_{0.889}^{0.983}$|	0.73 ± 0.33
Babi	DCB	23°36΄36″	75°39΄0″	71°34΄	59,399	2	|$1.045_{0.986}^{1.078}$|	1.02 ± 0.33
Bastet	SC	−2°3΄36″	20°27΄36″	49°51΄	77,026	8	|$1.007_{0.990}^{1.024}$|	0.92 ± 0.39
Bes	–	−54°43΄48″	115°19΄12″	93°54΄	120,370	9	|$0.983_{0.958}^{1.000}$|	−
Hapi	NCS	18°42΄36″	25°48΄0″	59°13΄	45,002	2	|$1.062_{0.993}^{1.088}$|	0.92 ± 0.33
Hathor	SC	25°4΄12″	25°21΄36″	48°18΄	45,184	8	|$1.019_{1.013}^{1.041}$|	0.98 ± 0.39
Hatmehit	D	3°17΄60″	−10°58΄48″	22°57΄	22,911	1	|$0.946_{0.887}^{0.989}$|	1.00 ± 0.39
Imhotep	NCS	−14°40΄48″	140°9΄36″	76°7΄	302,735	10	|$1.013_{0.998}^{1.031}$|	1.02 ± 0.27
Khepry	SC	−15°39΄36″	78°22΄12″	35°27΄	52,261	9	|$1.019_{0.998}^{1.035}$|	1.00 ± 0.31
Khonsu	–	−25°18΄36″	−153°33΄0″	79°45΄	57,202	1	|$1.013_{0.940}^{1.082}$|	−
Ma’at	NCD	32°16΄12″	9°24΄0″	37°54΄	57,365	2	|$1.003_{0.959}^{1.035}$|	1.00 ± 0.33
Neith	–	−27°18΄0″	−47°26΄24″	46°57΄	18,695	1	|$0.956_{0.886}^{1.012}$|	−

Similarly to this paper, Filacchione et al. (2016) have estimated the w for several morphological regions visible during the early phases of the mission. However, their absolute values are not comparable to ours as they lack observation α < 40°. Hence we concentrate on comparisons to their relative albedo, i.e. w/w₅₅₀ in Table 3. At first glance, all regions look shifted to darker ratios with respect to r_F/r_Hapke. Hapi and Hathor, two unambiguously bright regions, are slightly darker than average in Filacchione et al. (2016). We then rely on possible photometric correction mismatches as the explanation for such disparity in the albedo ratios.

6 DISCUSSION

6.1 ROI

When the phase curves of ROI and all-data are put together we observe an interesting example of the OE becoming progressively narrower with an increase in single-scattering albedo. In Fig. 7, the linear–exponential fits of each ROI and the all-data solution are compared and their exponential and linear terms are shown separately. We observe no large difference between the phase curves of the all-data and the Circular Mesa. ROI 1 is obviously integrated with the common albedo behaviour of the whole observed nucleus. With respect to the blue veins and the bright spots, it becomes clear that an incremental displacement of the phase curves to higher albedos and a sharpening of the exponential term is taking place. The HWHM of the blue veins is 12 per cent narrower than the whole nucleus, whilst for bright spots it is a further 37 per cent narrower. Their amplitudes (I) are 1.86 and 1.67, respectively, much lower than the all-data I of 2.86. Due to the reduction of the exponential term in describing both regions of interest, the linear term takes over, causing b and k to increase. The same parametric interplay is perceived among w, g_sca and h_SH for the Hapke modelling. Bright spots show less back-scattering and higher single-scattering albedo and sharper SHOE-H. Blue veins show intermediate values.

Figure 7.

Comparison of phase-curve morphologies. (a) Linear–exponential curves. (b) Exponential term, normalized to unity. (c) Linear term. (d) Linear term normalized to unity. The envelopes are the propagated uncertainties from the residuals and the Lommel–Seeliger correction (Appendix C). From the residuals we compute the offset and the standard deviation for bins of 2° of the phase angle. The deviation is larger for the measurements taken at 12° < α < 18° due to the oblique emergence angle.

The sharpening is therefore consistent with the incremental increase of some few per cent of ices mixed with the very dark component at the uppermost surface (Fornasier et al. 2016; Barucci et al. 2016). According to SHOE-H equation (Hapke 2012), h_s decreases as the packing factor decreases (i.e. superficial porosity increases). Therefore, if we consider that the ice grains permeate the dark grains and only widen the average inter-grain distance, h_s will decrease as ice abundance (or equivalent particle size) increases. The modification of superficial porosity by transparent water ice can also be interpreted as a reorganization in the deposition of dark grains in the optically active layer.

Concurrently, we cannot ignore some partial contribution of CBOE at very small phase angles for the ROI 2 and 3 due to the water-ice particles acting as internal scatterers. Mixtures of high- and low-albedo components have been proved to amplify second-order scattering and to weaken shadow-hiding (Shkuratov & Ovcharenko 1998; Nelson et al. 2004; Zubko et al. 2008). In this scenario, we have the average dark cometary surface being enriched by sublimating subsurface ice grains that increase the multiple scattering and CBOE effects, causing the slight OE sharpening we observe. However, the spatial resolution and average enrichment rate are not sufficient to further change the phase curve and produce the sharp spike at very small phase angles related to the CBOE as observed on the rings of Saturn, for example (French et al. 2007; Déau et al. 2013a).

The CBOE peak may only be indisputably verified if a single fresh bright feature is measured at very small phase angles at high spatial resolution. In this way, the bright spot signal will not be further diluted due to areal mixing inside the pixels. Nonetheless, some large albedos achieved by a few bright spots at α < 3° (see Fig. 6) may already point to the presence of this mechanism.

6.2 Comparison to laboratory measurements

Jost et al. (2017) used the PHIRE-2 radiogoniometer to measure the scattering curve before and after sublimation of inter/intra-mixture samples of 0.7 wt per cent 67 ± 31 μm ice particles, 0.2 wt per cent carbon black and 0.1 wt per cent tholins at 750 nm under 5° < α < 180°. The samples were fabricated on the SCITEAS simulation chamber (Pommerol et al. 2011, 2015a,b). The inter- and intra- mixtures differ in the structure of the sample at the scale of the individual grains. The mineral and water-ice grains are intimately mixed but the grains are individual units in the inter-mixture. In the intra-mixture preparations the mineral grains are contained as inclusions within water-ice particles. Jost et al. (2017) also measured the radiometry of a pure iceless mixture of tholins (33 per cent) and carbon black (66 per cent). Details of the preparation methods and the physicochemical properties of the mixtures are found in Poch et al. (2016a).

In this work, we have determined the phase curve of ice-rich features, which now provides a new comparison to the laboratory measurements. To verify any match between laboratory samples and cometary data, we decided to normalize all phase curves to α = 5° due to the spectral mismatch (Feller et al. 2016). We only selected measurements at i ≈ 0° and we concentrated our comparison on α < 16°, hence diminishing any possible contribution from macro-roughness due to large casting shadows. Then we traced a line from α = 5° to 16° to the Lommel–Seeliger-corrected A_eq of both cometary data and analogues. Fig. 8 presents the phase slopes of artificial cometary analogues along with those derived for all-data, blue veins and bright spots. In the figure, we have an explicit nonlinear trend for the flattening of the phase slope as the albedo increases. Higher albedo is a direct analogue for water-ice abundance on cometary surfaces. The desiccated tholin+carbon sample is darker than the nucleus, but tholin+carbon has a slight flatter phase slope. Bright spots, on the other hand, have the flattest phase slopes of our cometary data, approaching those measured for sublimated intra-mixtures. Sublimated intra-mixtures are almost completely dried out (Poch et al. 2016a,b); thus their closeness to the bright spots corroborates the very small water-ice abundances found by Barucci et al. (2016). Blue veins are slightly brighter than the all-data of the nucleus, but their phase slopes do not greatly differ. Off the small-phase-angle observations, blue veins are only distinguishable from the rest of the surface through spectrophotometry, which accentuates their peculiar phase curves.

Figure 8.

Linear phase slope of artificial cometary analogues (circles) along with those derived for all-data, blue veins and bright spots (triangles) with respect to their respective A_eq at α = 5°. Cometary analogues were measured at 750 nm; therefore they had to be spectrally recalibrated to 618 nm. The spectra are available in Poch et al. (2016b). The logarithm of the phase slopes is preferred for better visualization of the trend.

6.3 Optical layer structure

We compare our results to those obtained from instruments on-board the Philae lander. MUPUS (Spohn et al. 2015) and SESAME-PP (Lethuillier et al. 2016), respectively, estimated the thermal inertia (85 ± 35 J m⁻² K⁻¹ s^−1/2) and the electrical constant (2.45 ± 0.2) for a few centimetres to first metre layer of the Abydos site (Bibring et al. 2015). Both results point to a sintered microporous mantle covering the local area and possibly the whole comet to some extent. Near-surface porosity, constrained from about 30 per cent to 65 per cent for the MUPUS instrument and less than 50 per cent for SESAME-PP, may seem inconsistent with the photometric estimation of ∼80 per cent from the Hapke IMSA model (see Table 2); however, we must emphasize that the scales of depth are very different. For this, let us consider a rough estimation of mean radiative penetration depth D = λ/4πk for a semi-infinite continuum medium (Feynman 1964). In the hypothesis of extinction index k for amorphous hydrated carbon (k = 0.5) or tholin (k = 10⁻³) at 612 nm (Ciarniello et al. 2011) as analogues of optical indexes of the comet carbonaceous composition, we find that the optically active layer is constrained from ∼70 nm to ∼ 50  μm at most. In fact, this depth is far from the scales studied by the Philae instruments. The measured superficial porosity is then strictly related to irregularities in grains (Fray et al. 2016; Hilchenbach et al. 2016) and very top inter-grain spatial deposition (the ‘fairy castle’ structures of Hapke 1993). Simulations on fractal dust-aggregate properties and dynamics have proved such a low packing density to be possible (Levasseur-Regourd, Zolensky & Lasue 2008; Lasue et al. 2009).

6.4 Comparison to pre-perihelion results

The first all-data Hapke analysis of 67P/Churyumov–Gerasimenko using OSIRIS images was done by Fornasier et al. (2015) based on the pre-perihelion images up to 2014 August within 1.3° < α < 54°. Table 4 reproduces the Hapke parameters from their all-data Hapke (2012) analysis. It is noteworthy that the Fornasier et al. modelling was undertaken at 649 nm, whilst we only have data at 612.8 nm. By verifying the spectral slopes measured in their paper, the Imhotep region and its vicinities present a slope of about 12.5 per cent/(100 nm) for α < 10°. As such, we recalibrated all normal and single-scattering albedos by 4.5 per cent. The new values are listed in Table 4 under w₆₁₃ and ρ_ν,613.

Table 4.

Hapke IMSA parameters previously obtained for the nucleus of 67P/Churyumov–Gerasimenko.

	w649	w613	B0	hs	gsca	Porosity	|$\bar{\theta }$|	ρν,649	ρν,613
Fornasier et al. (2015)	0.034	0.0325	2.25	0.061	−0.42	0.87	28	6.7 per cent	6.4 per cent
Feller et al. (2016)	0.038	0.0355	2.56	0.067	−0.37	0.86	15.6	6.8 per cent	6.5 per cent
‘Bright spot’	0.067	0.063	2.07	0.103	−0.26	0.81	[15]	7.7 per cent	7.4 per cent
‘Sombre boulder’	0.029	0.027	2.27	0.108	−0.41	0.82	[15]	6.4 per cent	6.1 per cent

The aim of such a comparison is to check for any substantial darkening or brightening on the surface of the nucleus after the perihelion passage. At first glance of the Hapke parameters, we observed that w_all-data and h_{s, all-data} are different by −0.055 (−17 per cent) and 0.02 (+33 per cent) compared to w_{613,Fornasier} and h_s,Fornasier. This dichotomy is higher than the associated errors and has as a possible explanation the presence of a bright water-ice-rich region in the field-of-view. Bright regions elevate the average apparent albedo, elevating w and broadening h_s. Our data are mostly constrained to Imhotep and its vicinities whereas all images analysed by Fornasier et al. have the bright Hapi, Hathor and Seth at α < 15° (see fig. 5 in Fornasier et al.) in the field-of-view. Seth is absent in our data, whilst Hapi and Hathor are partially present in just two images. This is due to the fact that the Southern hemisphere became more visible for a few months before and after the perihelion passage, unlike in 2014, when the Northern hemisphere was visible and better insolated. Therefore, the all-data Hapke parameters obtained in this paper are representative of a desiccated surface, mainly Imhotep and its surroundings.

To be able to check for actual albedo alteration, we chose a single image taken at |$\bar{\alpha }=1.37^{\circ }$| on 2014 July 29 UT 00:45:31 (r = 3.65 AU) that exhibits most of Imhotep, Ash and Aten to compare to the 2016 April 9 UT 23:59:32 image (Fig. 3). Both images were photometrically corrected using the Hapke model and then normalized to e = 30° and i = 30° (Appendix B). The colour gradient represents the absolute equigonal albedo at |$\bar{\alpha }=1.37^{\circ }$|⁠. The 2016 April 9 UT 23:59:32 image was projected onto the 2014 July 29 observational configuration for better visualization. Both are shown in Figs 9(a) and (b). When comparing them, we observed a reduction of a central bright feature on Ash along with a mild darkening (approx. −2 per cent) in respect to 2014. Imhotep otherwise seems to have experienced no outspread darkening (<0.5 per cent) over its dusty deposit. The darkening is possibly connected to seasonal depletion of water on buried layers or redeposition of dust grains through cometary activity.

Figure 9.

Checking for post- and pre-perihelion albedo alteration. (a) Hapke photometrically corrected image at 2014 July 29 UT 00:45:31 (⁠|$\bar{\alpha }=1.37^{\circ }$|⁠) using parameters from Fornasier et al. (2015) (Table 4). (b) Hapke photometrically corrected image at 2016 April 9 UT 23:59:32 projected onto the 2014 July 29 UT 00:45:31 observational configuration. (a) and (b) are normalized into an equigonal albedo corresponding to e = 30° and i = 30° (Appendix B) and to λ = 613.8 nm. (c) Hapke phase curves for 2016 April 10 (in blue) and 2015 February 14 (in red, from Table 4) of the Imhotep–Ash flyby area studied by Feller et al. (2016). Data from this paper are underplotted in grey-scale.

Finally, to check whether we identify any seasonal albedo variation at the metre-scale, we compared the phase curve we retrieved for the same area visited by the 2015 February 14 flyby (Feller et al. 2016). The whole flyby area is about 200 metres and is located on the Imhotep–Ash boundary. Similarly to the comparison to Fornasier et al. (2015), we had to spectrally recalibrate the Hapke parameter of Feller et al. (2016) using a slope of 17.7 per cent/(100 nm) at 649 nm. It corresponded to a spectral shift of 6.4 per cent (see the values in Table 4). The area shows albedo variations <4 per cent, in good agreement with Feller et al. (<3.5 ± 2.1 per cent). To compose its phase curve, we retrieved all shape-model facets under latitude = {4°21^΄, 17°1^΄} and longitude = {−170°17^΄, 178°58^΄} from the 2016 April and February 10 sequences. We then fitted the data with the Hapke IMSA model (Appendix B) and obtained the following parameters: w = 0.0255 ± 0.002, g_sca = −0.41 ± 0.02, B₀ = 2.62 ± 0.3, h_s = 0.095 ± 0.01 and |$\bar{\theta }=39^{\circ }\pm 5^{\circ }$|⁠. The parameters are indeed very different from those obtained in 2015 February, but strikingly similar to the ‘sombre boulder’ found in the area (Table 4). As shown in Fig. 9(c), where the phase curves associated with each set of Hapke parameters are plotted together, the area seems to have become darker and its SHOE-H, broader. At 618 nm, the normal albedo was expected to reach 6.4 per cent according to the previous Hapke parameters, but is no larger than 6.0 per cent instead. The area seems to have been affected by the same process that mildly darkened Ash. Conversely, the difference of 0.5 in albedo indicates a darkening of about −7 per cent, which is far larger than that measured for Ash. This could represent a particular intensification of the process at a local scale, where airfall dust may accumulate at certain topological features.

7 CONCLUSION

The OSIRIS camera on-board Rosetta has obtained images of the nucleus at very small phase angles for the third time. Acquired on the night run of 2016 April 9–10, the observations show the OE on the filter WAC F18 sweeping across Imhotep with phase angles ranging from 0° to about 17°. To better constrain the phase-curve morphology, we have added the 2016 February 10 sequence, which also displays Imhotep at full extent. From our photometric analysis of all-data, a Circular Mesa (ROI 1), a blue vein (ROI 2) and an assemblage of bright spots (ROI 3), we measured a stepwise SHOE-H narrowing from the reddest (ROI 1, ∼18 per cent/100 nm) to the bluest (ROI 3, <9 per cent/100 nm, Oklay et al. 2016b; Deshapriya et al. 2016; Fornasier et al. 2017) colours. The HWHM goes from 7.27° to 4.27° as the single-scattering albedo increases from 6.14 per cent to 7 per cent, measured by our linear–exponential fitting. Therefore we interpret this behaviour as an outcome of the increase in the abundance of ices inter-mixed with dark grains on the blue veins (ROI 2) and bright spots (ROI 3). The trend of linear phase slopes for cometary analogue mixtures and our cometary data shows a smooth transition from desiccated to water-rich cometary surfaces. Here, we observe bright spots having similar phase slopes to sublimated intra-mixtures (Poch et al. 2016a,b), which corroborates the very small water-ice abundance found by Barucci et al. (2016).

Once the photometric parameters were retrieved, we undertook a photometric correction of the full 2016 April 10 set. We measured the albedo ratio for all cometary regions observed therein. In particular, we found Hapi to be the brightest, up to +8.8 per cent higher than average, while Hatmehit is shown as the darkest, with −8.9 per cent at minimum. The far better insolation of the Northern hemisphere, where the bright regions such as Hapi, Hathor and Seth are found, is sufficient to explain mismatches on the Hapke parameters from pre- and post-perihelion epochs (Fornasier et al. 2015). Next, we compared one of our photometrically corrected images to a pre-perihelion image taken at the smallest phase angle (⁠|$\bar{\alpha }=1.37^{\circ }$|⁠) on 2014 July 29 UT 00:45:31 (Fornasier et al. 2015) to check for albedo alteration. We found only mild darkening (−2 per cent) on the Ash region. On the other hand, for the phase curve of the 2015 February 14 flyby area (Feller et al. 2016) we observe a darkening of up to −7 per cent in normal albedo. When we put our phase curve next to the one from Feller et al. (2016), it is remarkable that the OE became flatter with respect to 2015. Darkening may have locally affected some certain areas where airfalling of dust was favourable, whilst altogether the global albedo remained unchanged.

Acknowledgements

We thank Dr Karri Muinonen for his review and thoughtful comments. The authors thank all the developers and keepers of Python libraries such as Scikit-image, Pandas, Astropy and Numpy. We also thank Matej Tyc and Christoph Gohlke for implementing and making available an image registration code based on the discrete Fourier transform.

We acknowledge the financial support by Île de France (DIM ACAV). OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung, Göttingen, Germany; CISAS–University of Padova, Italy, the Laboratoire d’ Astrophysique de Marseille, France; the Instituto de Astrofísica de Andalucia, CSIC, Granada, Spain; the Research and Scientific Support Department of the European Space Agency, Noordwijk, The Netherlands; the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain; the Universidad Politéchnica de Madrid, Spain; the Department of Physics and Astronomy of Uppsala University, Sweden; and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB), and the ESA Technical Directorate is gratefully acknowledged.

1

We have attached H to the shadow-hiding acronym to differentiate the actual phenomenon from its mathematical representation given by B_SH in the Hapke model (see Appendix B).

2

http://www.cosmos.esa.int/web/spice/spice-for-rosetta

3

IMREG_DFT package, http://imreg-dft.readthedocs.io/en/latest/

APPENDIX A: LINEAR–EXPONENTIAL FITTING

A simple linear–exponential formula is applied to obtain morphological parameters for photometric phase curves:

\begin{equation} r_\textrm{l--e}(\alpha ,\lambda )=A\cdot \exp (-\alpha /d) + b-k\cdot \alpha , \end{equation}

(A1)

where A is the amplitude of the opposition surge, d is the width of the opposition surge, b is the phase-curve albedo without opposition surge and k is the angular coefficient of the linear part. Accordingly, the interpretation of these parameters would be model-free and easily comparable to those available to other surfaces or small bodies (Kaasalainen et al. 2003). The associated HWHM is then calculated as (Rosenbush et al. 2005):

\begin{equation} \mathrm{HWHM}=d/1.45. \end{equation}

(A2)

For correcting the brightness trend from limb to terminator and also the topographic–photometric effects, we applied the Lommel–Seeliger law (Fairbairn 2005). This law comes directly from the radiative transfer equation and is adequate as disc function for dark, single-scattering surfaces, such as the comet nucleus:

\begin{equation} D_\mathrm{LS}(i,e)=2\frac{\mu _{0}}{\mu +\mu _{0}}. \end{equation}

(A3)

The topographic–photometric corrected radiance factor at a mirror point fixed at |$\bar{\alpha }$|⁠, the so-called equigonal albedo (Shkuratov et al. 2011), is then simply calculated through

\begin{equation} A_\mathrm{eq}(\bar{\alpha },\lambda ^{\prime })=\frac{r_{F(\alpha ,i,e,\lambda )}}{r_\textrm{l--e}(\alpha ,\lambda )\cdot D_\mathrm{LS}}\cdot r_\textrm{l--e}(\bar{\alpha },\lambda ^{\prime }), \end{equation}

(A4)

where μ = cos (e), μ₀ = cos (i) and the equigonal albedo is equivalent to |$r_{F}(\bar{\alpha },\bar{\alpha }/2,\bar{\alpha }/2,\lambda )$|⁠.

In an attempt to interpret the parameters of very dark surfaces as in the case of 67P, A and d are mostly connected to the SHOE amplitude and width and thus properties such as packing factor, micro-roughness, particle irregularities and opacity (Hapke 1993), whereas b and k are partially related to the single-particle phase function and also roughness in different scales (Kaasalainen 2003; Helfenstein & Shepard 2011).

We then relied on the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS, Broyden 1970; Zhu et al. 1997) in the basin-hopping method to minimize the density-weighted residual root mean square (RMS) between the facet radiance factor and equation (1). The basin-hopping method has been widely used to solve complex molecular structures (e.g. Wales & Doye 1997; Verma et al. 2006; Prentiss, Wales & Wolynes 2008) and works by scanning several local minima though random controlled displacements in the parametric space until the global minimum is found, i.e. the real solution. The BFGS is used for local minimization inside the basin-hopping method. To secure the complete independence of the parameters and consolidate the best solution as the global one, we selected 20 initial conditions of lowest RMS over 1000 randomly simulated sets of parameters. The parametric uncertainties were computed through the RMS-weighted standard deviation of all solutions obtained.

Déau et al. (2009) and Kaasalainen et al. (2003) have remarked that the linear–exponential formalism appears to depend on phase-angle coverage. As a matter of fact, to properly retrieve stable parameters the data must be sampled at three critical points: the surge rising near zero phase angle, the linear-to-exponential curve concave, and some scatter under the linear part (15° ≲ α ≲ 60°). We satisfactorily span these three points with the addition of the 2016 February 10 sequence into our data. All studied phase curves, from all-data to ROI, are sampled evenly, i.e. phase-angle coverage is consistent throughout our analysis.

APPENDIX B: POROSITY-DEPENDENT HAPKE IMSA MODEL

Hapke (2002, 2008, 2012)) has proposed in the last decade an improved treatment of the role of filling factor and of multiple scattering in the IMSA (isotropic multiple scattering approximation) model. This version of the model (Hasselmann et al. 2016) has previously been applied to the nucleus (Fornasier et al. 2015; Feller et al. 2016), retrieving parameters compatible to previous results on cometary surfaces. Since the data are incomplete for α > 17°, free parameters such as average macroscopic roughness slope (⁠|$\bar{\theta }$|⁠) and asymmetry factor (g_sca) may have their values better constrained when images at α ≫ 80° are analysed (Schmidt & Fernando 2015). Other variables, i.e. shadow-hiding amplitude B₀, shadow-hiding width h_s and porosity factor K, are otherwise dependent on very small phase-angle measurements. The bi-directional reflectance r_Hapke, disregarding Akkermans et al. (1988) formalism for the CBOE, is described by the following equation:

\begin{eqnarray} r_{\mathrm{Hapke}(\alpha ,\mu _{0},\mu ,\lambda )}&=&K_{(h_{s})}\frac{w_{\lambda }}{4}\left(\frac{\mu _{0e}}{\mu _{e}+\mu _{0e}}\right)[(1+B_{\rm SH(\alpha )})p_{(\alpha )}\nonumber\\ &&+\,(H_{(\mu _{0}/K,w_{\lambda })}H_{(\mu /K,w_{\lambda })}-1)]S_{(\mu _{0},\mu ,\alpha )}, \end{eqnarray}

(B1)

where μ_0e and μ_e are the effective cosines of the incidence and emergence angles, involving the topographic correction of the facet by a macro-roughness shadowing function S. w_λ is the single-scattering albedo, B_SH is the shadow-hiding OE term as described in Hapke (1993), p is the single-lobe Heyney–Greenstein particle phase function (SPPF) and H is the isotropic multiple-scattering function, analytically described by the second-order approximate Ambartsumian–Chandrasekhar function (Hapke 2002). A further qualitative and mathematical description of this applied version of the model and functions therein can be found in Helfenstein & Shepard (2011) or the Hapke (2012) book. The minimization technique and uncertainty estimation are the same as those adopted for the linear–exponential formalism in Appendix A.

Once the Hapke parameters are obtained, the topographic–photometric correction is then carried out as follows:

\begin{equation} A_\mathrm{eq}(\bar{\alpha },\lambda ^{\prime })=\frac{r_{F(\alpha ,i,e,\lambda )}}{r_{\mathrm{Hapke}(\alpha ,\mu _{0},\mu ,\lambda )}}\cdot r_{\mathrm{Hapke}(\bar{\alpha },...,\lambda ^{\prime })}, \end{equation}

(B2)

where A_eq is the equigonal albedo. The standard procedure is to normalize A_eq into the e = 30° and i = 30° configuration.

APPENDIX C: UNCERTAINTIES ASSOCIATED WITH THE LOMMEL–SEELIGER CORRECTION

The Lommel–Seeliger law has been widely used for disc correction in the scope of the spectrophotometric analysis of OSIRIS images (Fornasier et al. 2015; La Forgia et al. 2015; Lucchetti et al. 2016; Oklay et al. 2016b; Barucci et al. 2016; Pajola et al. 2016; Deshapriya et al. 2016; Feller et al. 2016; Oklay et al. 2016a; Fornasier et al. 2016). Therefore, we use this appendix to briefly discuss the role of multiangular uncertainties in the estimation of A_eq. We remind the reader that the (i, e, α) angles are product of the facets’ normal vectors, thus their imprecisions are associated with the shape-model imprecision itself. From the Lommel–Seeliger law (equation 3), we may apply the error propagation formula to derive the relative error |$\eta _{A_\mathrm{eq}}$|⁠:

\begin{equation} \eta _{A_\mathrm{eq}}^{2}=4\left[\left(\frac{\eta _{r_{F}}}{D_\mathrm{LS}}\right)^{2}+\left(\frac{\upsilon _{\mu }}{\mu _{0}}\right)^{2}+\left(\frac{\mu }{\mu _{0}}\upsilon _{\mu _{0}}\right)^{2}\right], \end{equation}

(C1)

where η represents the relative errors associated with the subscript quantities, whilst υ represents the absolute errors. |$\eta _{r_{F}}$| is already discussed in section 2. υ_μ and |$\upsilon _{\mu _{0}}$| are otherwise unknown and would depend on a deep comparison of the shape model to the high-phase-angle images. Based on equation (6), we report (i, e) maps for |$\upsilon _{\mu ,\mu _{0}}=\pm \left\lbrace 0.1^{\circ },0.5^{\circ },1^{\circ },5^{\circ }\right\rbrace$|⁠, as shown in Fig. C1. We observe two regimes in the (i, e) gradients: when |$\upsilon _{\mu ,\mu _{0}}$| weighs less than |$\eta _{r_{F}}$| (±0.1° and ±0.5° maps), the uncertainties are mostly independent of the incidence angle for i < 50°; and when |$\upsilon _{\mu ,\mu _{0}}$| becomes the main error source, the gradient profile is radial. In either case, all disc-corrected measurements over i > 75° are largely unreliable. But the same divergence is not reproduced over the emergence angles. |$\eta _{A_\mathrm{eq}}$| are kept under about 10 per cent when i < 60° for |$\upsilon _{\mu ,\mu _{0}}\le \pm 1^{\circ }$|⁠. Through visual inspection of the images with and without disc correction, we do not detect any large disparity among the albedo features; thus we estimate |$\upsilon _{\mu ,\mu _{0}}\le \pm 1^{\circ }$| (⁠|$\bar{\delta }_{A_\mathrm{eq}}\lesssim 5{\rm {per\,cent}}$|⁠) for the most part of the normal vectors.

Figure C1.

(i, e) maps represent the A_eq uncertainty gradients due to Lommel–Seeliger disc correction. The maps are, respectively, provided for |$\upsilon _{\mu ,\mu _{0}}=\pm \left\lbrace 0.1^{\circ },0.5^{\circ },1^{\circ },5^{\circ }\right\rbrace$|⁠. The uncertainties are far more sensitive to incidence angle than emergence angle.

Regarding the uncertainties associated with the phase angle, it is solely dependent on altitude and image registration error. Through the propagation formula, we have:

\begin{equation} \upsilon _{\alpha }^{2}\approx \arccos \left(1-\frac{\delta _{XY}^{2}}{2D^{2}}\right), \end{equation}

(C2)

where δ_XY is the absolute translation error from image registration and D is the sub-spacecraft altitude. Therefore, |$\upsilon _{\alpha }^{2}\approx 0.0055^{\circ }$| when considering a quite high |$\bar{\delta }_{XY}=\pm 3$|px (OSIRIS–OASIS image subtraction gives |$\bar{\delta }_{XY}<3$|px off the image borders), which is much smaller than the figure ticks.