© The Authors 2024

1 Introduction

Phobos and Deimos, the moons of Mars, are quite peculiar in relation to other moons in the Solar System. They are very small compared to Mars (~15–25 km in size), and they have heavily cratered surfaces and a relatively low density. Phobos, in particular, displays a complex geomorphology, characterized by some peculiar 1 km long grooves. Several of its craters show various states of degradation. Notably, the Stickney crater has a diameter of up to 9.4 km (Thomas et al. 1999). Phobos is covered by a layer of regolith 5–100 m thick (Basilevsky et al. 2014). Its surface is dark, as it is for Deimos, with a low geometric albedo (~7%) in the υ filter (Zellner & Capen 1974). Its bulk density is about half of Mars (1.85±0.07 g cm⁻³, Pieters et al. 2014).

The surface composition is heterogeneous, and two distinct units, blue and red, have been identified since the first spatially resolved spectra acquired by the Phobos 2 mission (Murchie & Erard 1996; Bibring et al. 1992). The blue unit shows a moderate slope spectrum and higher reflectance, while the red unit is darker and has a steeper slope, especially in the near-infrared region. The red unit dominates the majority of the surface of Phobos, while the blue one is mostly located around and partially inside the Stickney crater (it was initially considered a fresh ejecta deposit). Further high-resolution images acquired with the HIRISE instrument on board the Mars Reconnaissance Orbiter (MRO) indicated a more complex mixing in composition, with spots of redder material within the blue unit areas (Basilevsky et al. 2014).

The origin of the satellites of Mars is debated, and two theories have been proposed. One theory is that Phobos and Deimos are asteroids captured by Mars (Hansen 2018), a hypothesis reinforced by the spectral similarity between the two satellites and primitive D-type asteroids, which are found mainly in the outer asteroid belt and among Jupiter’s Trojans. However, dynamical simulations of the possible gravitational capture scenario fail to explain their quasi-equatorial orbits (Burns 1992). The spectral analogy with D-type asteroids suggests that the Martian satellites, if they are indeed captured asteroids, may be rich in organic material, carbon and hydrous silicates, and may possibly contain water ice in their interiors (Rivkin et al. 2002). The second theory claims that the satellites formed in situ from a debris disk produced by the collision of a large object with Mars (Craddock 2011; Rosenblatt & Charnoz 2012). This scenario explains the orbits of these moons, but it does not explain their carbonaceous-like compositions (Nakamura et al. 2021).

The spectra of Phobos are usually featureless, except for a 0.65 µm absorption feature observed in the red unit and a sharp 2.8 µm feature observed in both units (Fraeman et al. 2014). The red unit feature centered at 0.65 µm was observed by a ground-based telescope (Murchie et al. 2015) and from different mission data (Simonelli et al. 1998; Fraeman et al. 2014; Pajola et al. 2013), and it has been attributed to nontronite, a desiccated iron bearing phyllosilicate (Murchie et al. 2008), or to space weathered anhydrous silicates (Clark et al. 2012). The 2.8 µm feature is stronger in the red unit spectra and attributed to structural hydroxyl embedded in mineral lattices (Rivkin et al. 2002; Fraeman et al. 2014).

The debate on the origin of the Martian satellites should soon be resolved thanks to samples that will be returned by the JAXA Martian Moons eXploration (MMX) mission, which is devoted to the exploration of the Mars system. The MMX mission will be launched in 2026, be inserted into Mars orbit in 2027, and investigate the Martian system for three years, focusing mainly on Phobos, the key target of the mission. The main goals of MMX are to return samples of Phobos and, through both a detailed in situ investigation of the Mars satellites and further laboratory studies of Phobos samples on Earth, to clarify the origin of the Mars satellites and the process of planet formation in the Solar System.

In this work, we perform the very first analysis of Phobos photometric properties using data from the High Resolution Stereo Camera (HRSC) on board the Mars Express spacecraft. This work will support the MMX mission planners in the optimization of the trajectories and science operations, especially in defining the exposure time needed to achieve the expected signal-to-noise ratio of a given Phobos observation.

Photometric investigations of Phobos have been carried out in the past by Klaasen & Duxbury (1979), Pang et al. (1983), Avanesov et al. (1991), Simonelli et al. (1998), Cantor et al. (1999), and Thomas et al. (1999) and more recently by Pajola et al. (2012) and Fraeman et al. (2012). The HRSC data offer a unique opportunity to study in detail the Phobos opposition surge, that is, the non-linear increase of the reflectance for small phase angle values.

2 Observations and data reduction

We investigated Phobos’ photometric properties using the Mars Express HRSC data. The Mars Express mission was launched in 2003 by the European Space Agency, and it is principally devoted to the study of the Martian surface and atmosphere, especially to the observation of traces of water on the planet’s surface.

The HRSC camera operates as a push-broom scanning instrument with an integrated stereo capability. It is used to observe the Martian surface and atmosphere as well as the planet’s two moons, Phobos and Deimos. This advanced camera system comprises nine CCD line detectors that capture highly detailed, simultaneous high-resolution stereo image swaths, and it has the ability to incorporate multiple colors and phase angles (Neukum & Jaumann 2004; Jaumann et al. 2007). Among the nine line sensors, four are equipped with spectral filters, enabling the acquisition of images in red, green, blue, and infrared wavelengths. The remaining five panchromatic sensors (nadir, stereo 1 + 2, photometry 1 + 2) serve to determine the photometric properties of the surfaces of the celestial bodies and enable the creation of digital terrain models (DTMs) (Neukum & Jaumann 2004; Gwinner et al. 2009), whereas DTMs are only produced from observations of the Martian surface. The HRSC data calibration involves radiometric calibration for all acquired images and geometric correction for Mars surface observations.

The additional super resolution channel (SRC) of HRSC operates in panchromatic mode, utilizing a single CCD detector to acquire images with a significantly higher pixel resolution compared to the other color channels of the camera. The SRC captures images with a pixel size of 2.4 meters at an altitude of 250 kilometers (Jaumann et al. 2007) and was initially intended to allow for the detailed examination of geological features, surface textures, and small-scale structures on Mars. With higher resolution images available from other Mars orbiter missions, the main focus of the SRC turned to the observation of the Martian moons Phobos and Deimos as well as astrometric observations. A typical calibration process for SRC could not be realized due to time and budget constraints during the mission preparation phase. In addition to the usual functional tests during the development and implementation process of the electronics, the only radiometric calibration measurements have been flat-field determinations by means of an integrating sphere and dark current measurements (DSNU). The calibration process of the SRC image data consists of four steps: (1) calculation of the average of the dark current pixels at the border of the CCD to get the dark signal uniformity (DSU) level, (2) subtraction of the product DSU×DSNU from the image, (3) division of the image by the flat field, and (4) the flagging of hot pixels. This process permits the elimination of hot pixels generated on the CCD sensor by solar flares during the Mars Express cruise phase (Jaumann et al. 2007).

Phobos photometry observations are acquired whenever the observation opportunities meet favorable requirements in terms of spatial resolution and small phase angle coverage, that is, when Phobos is observable at relatively small phase angles (<10°) and when the Mars Express-target distance is less than 10 000 km. Depending on the distance to the moon, either SRC observations with inertial pointing (at higher distances) or HRSC color observations (blue, green, red, and infrared filters) with a spacecraft slew pointing mode (at lower distances ~6000 km) are planned. For very small phase angles (<1°), SRC observations in spot tracking mode are used, where the spacecraft keeps the moon constantly in the center of the field of view. The resulting pixel resolutions range between 2.4 meters and 400 meters, depending on the observation distance and camera system used (i.e., HRSC or SRC).

From the ESA Planetary Science Archive, we first selected images acquired with the absolutely calibrated blue, green, red, and IR HRSC filters. These data cover mostly the 10–100° phase range, with very few observations at small phase angles.

Therefore, to study the opposition effect, we also analyzed the SRC observations.

2.1 HRSC data analysis

For the HRSC data in the absolutely calibrated blue, green, red, and IR filters, the radiance factor (also named I/F) was simply determined using the correction factors provided in the header of each image, which consider the absolute calibration values for each filter and the varying heliocentric distance. The integrated flux was evaluated from aperture photometry divided by the projected Phobos surface, that is, considering the illuminated pixels and those in shadows derived for a given observation by the simulated images. To mimic the geometric conditions, namely, incidence (i), emission (e), phase angle (α), longitude (lon), and latitude (lat), for each pixel, we generated simulated images using the latest version of the Mars Express mission and instrument spice kernels (including the latest shape, size and rotational constants¹). For the shape model, we utilized the stereopho-togrammetric shape model² generated from Mars Express HRSC images (Willner et al. 2014). For the simulation, we used the UTC mean time, that is, the time corresponding to the start time of a given acquisition plus half of the exposure time.

However, given the non-optimal NAIF/SPICE kernels, the generation of synthetic images was hindered. In fact, for several observations the simulated images were not in the nominal field of view of the camera, and non-systematic shifts in the field of view were needed. To improve the simulations, we adopted some field of view corrections files provided by the HRSC team that take into account the non-linearity of the pixel positions, but even with these corrections, some manual vertical and/or horizontal adjustments were still needed for a number of images.

In addition to this problem, we also discarded some simulations because their quality was not satisfactory. In fact, images were acquired in push-broom mode, and in some cases, the simulated images looked elongated and/or did not match completely with the original Phobos shape observed in a given image (Fig. A.1).

We also had to apply image registration to improve the pixel matching between real and synthetic images. Given that HRSC image acquisitions are based on the push-broom mode, the instrument induces occasional distortions to the imaged cross section of the body. Moreover, the synthetic images were hindered by the NAIF/SPICE kernel’s imprecision, which prevented the synthetic imaging from actually accounting for all distortion effects seen in the original images. In a past analysis performed by our team on space mission data, typical image co-registration of disk-resolved images were performed, including corrections in offset, bi-dimensional rotation, and projective transformation (Fornasier et al. 2015; Feller et al. 2016; Hasselmann et al. 2017). However, the mitigation of HRSC distortions must be dealt with by using a different technique. We therefore used the latest implementation of the TV-L1 dense optical flow technique (Perez et al. 2013). It computes the full flow field of displacement vectors with respect to changing features for pairs of images obtained under limited context alteration, that is, an entire body appearing in two images conserves most of the features in both frames. This estimated flow vector field transforms the synthetic images into similar pixel distributions of the original images, thus setting a correspondence between every pixel in the synthetic image with those in the original image. In our implementation, we fed the TV-L1 optical flow algorithm, available in the Python scikit-image package (Van Der Walt et al. 2014), synthetic images where Phobos’ surface brightness was estimated using the Lommel-Seeliger law. The synthetic images were compared to original HRSC images that had been signal normalized in order to obtain the flow vector fields. Then, the transformed synthetic images were visually inspected in order to tune up the free parameters in the algorithm. To estimate the registration quality, we sorted the Universal Quality Image Index (UQI; Wang & Bovik 2002), which calculates the amount of similarity between two images. The UQI has been reported to perform better against noise as well as blurriness. When compared to visual inspection, we set UQI > 0.5 as a satisfactory condition for similarity. The final list of images on which we based our analysis is shown in Tables B.1–B.4, available at CDS.

2.2 Phobos opposition surge: SRC data analysis

The SRC observations cover the Phobos opposition surge during six orbits from 2019 to 2022 (orbit numbers: K076, K122, K623, L132, L155, and M 247). However, we ultimately did not use the K623 observations because they are partially contaminated by Mars straylight. In fact, during this observational session, Mars Express was on the dayside of Mars, in contrast to its position when the data were acquired during the other orbits.

These observations were optimized to investigate the opposition surge, covering the phase angle range 0.98–17.10° for K076, 0.4–7.0° for K122, 0.24–2.96° for L132, and 0.10–6.70° for M 247, with a number of individual observations per orbit between 40 (K076) and 100 (L132). The final list of SRC images used to investigate the opposition surge is shown in Table B.5, available at CDS.

Examples of different phase angle images acquired in orbits K076 and M 176 are presented in Fig. 1, which clearly displays the increase of brightness as the phase angle decreases.

Due to a lack of absolute calibrations, we trust only in relative photometry for SRC images. Data were corrected applying the following procedure:

(a) We first divided each image by the exposure time because the data are provided in digital units and not normalized by the exposure time. (b) We corrected the data by the heliocentric distance (dist_⊙) to the square, similar to what is done when generating the radiance factor, in order to consider the varying incoming solar flux. (c) We computed the disk-integrated photometry for each observation simply as the integral of the flux from Phobos, evaluated using aperture photometry, divided by the object projected surface. The projected surface was estimated from simulated images generated using the latest Mars Express mission and instrument spice kernels available on the ESA planetary science archive and using the highest resolution Phobos shape model (Ernst et al. 2023).

In contrast to the HRSC calibrated filter, the simulated images of the SRC framing camera match pretty well with the shape of the original ones. The co-registration was also quite straightforward, as only a simple translation of a few pixels was needed to match the real images.

Some of the observations cover extremely small phase angles that are smaller than the angular size of the Sun as seen from Phobos: $${\gamma }_{\text{sun}}=\arcsin \frac{R_{\odot }}{{\text{ dist }}_{\odot }},$$(1)

where R_⊙ is the radius of the Sun and dist_⊙ is the heliocentric distance of Phobos at a given observation. The angular size of the Sun is 0.16° for orbit M247 (dist_⊙ = 1.66 au) and 0.19° for orbit L132 (dist_⊙ = 1.39 au). For these two orbits, which cover the smallest phase angle, we considered in the photometric analysis only the observations acquired at α > γ_sun to avoid an underestimation of the opposition effect. In fact, it is well known that the brightness values flatten for α < γ_sun (Déau 2012).

The integrated flux, in arbitrary units, is presented in Fig 2. All the data display a similar phase function behavior. In fact they are superposed after the normalization process described in the following (Fig. 3). As the data are not absolutely calibrated, we attempted to estimate the I/F SRC data as follows, even if we stress that the I/F shown in Fig. 3 is not an absolutely calibrated value. We first normalized the data in count since the observations during the different orbits were not at the same flux intensity. For this, we interpolated each orbit’s data at a finer phase angle resolution, and then we normalized all the observations at α = 2°. Finally, we multiplied the SRC data by a factor in order to match the calibrated reflectance values of the green filter at α = 5°. This normalization factor is of course not perfect because the observations were acquired at different spectral ranges. In fact the SRC spectral range is larger than the green one. Moreover, the observed surface of Phobos is not the same. We point out that we utilized the SRC data to derive the width and amplitude of the opposition surge and not to derive Phobos’ absolute reflectance. The estimated SRC I/F is presented in Fig. 3.

Fig. 1 Images showing the flux increasing with a decreasing phase angle from images acquired during the orbit K076 (top panel) and M 247 (bottom panel) with SRC. The flux is in DU s−1; no other normalization was applied here (the heliocentric distance was the same for a given orbit).

Fig. 2 Observations of Phobos opposition effect from five individual SRC observations. Data have been normalized by the exposure time and corrected for heliocentric distance to the square. The insert shows very small phase angle data acquired for orbits L132 and M 247. Data acquired at a phase angle lower than the angular size of the Sun seen from Phobos distance were discarded in the modeling.

3 Disk-averaged photometry

The disk-averaged photometry represents the global analysis of Phobos surface photometric properties. To model the Phobos photometry, we used the Hapke formalism following the same methodology and equations presented in Fornasier et al. (2015). For the global photometry, we used the disk-integrated Hapke expression (Hapke 1993) with a single-term Henyey-Greenstein function and neglecting the coherent-backscattering mechanism (CBOE), which is expected to have a low contribution for dark surfaces (Shevchenko et al. 2012). This model has five-parameters: the single-scattering albedo (SSA) w_λ, the asymmetry factor 𝑔_λ, the average roughness angle $\overline{\theta }$, and two parameters defining the amplitude B₀ and the width h of the shadow-hiding opposition surge.

Since observations at small phase angles are available only with SRC, we used the estimated SRC radiance factor on the combined observations from five orbits to determine the Phobos opposition surge parameters. As the Hapke parameters are not easy to disentangle and the roughness parameter may be determined only with observation at α > 60°, we fixed it at 24°, the value determined from the modeling of the disk-resolved images obtained with the green filter (see Sect. 4). We used the Levenberg-Marquardt algorithm (Moré 1978) to find the best fit between the Hapke model and the observations, looking for the minimum inside the boundaries w_λ = {0.001, 0.3}, 𝑔_λ = {−1.0, 1.0}, B₀ = {0.2, 5}, and h = {0.0, 0.5}. The uncertainties in the parameters were calculated from the covariance matrix in the minimization algorithm, and they were derived from the square root of the diagonal elements of this matrix.

The results are reported in Table 1. Since we fixed the roughness value, we tested the stability of the results reported in Table 1 for $\overline{\theta }$ ranging from 21° to 27°, finding that for the lower and higher roughness limits, w and 𝑔 slightly increase and decrease, respectively, but with values still within the uncertainties reported in Table 1.

Since the SRC data are not absolutely calibrated, we did not report the SSA in Table 1 to avoid confusion. The estimated fitted value is w = 0.088 for the SRC data, higher than what is obtained for the green and red filters (Table 1) but lower than the value derived for the IR filter. Considering the very red spectrum of Phobos, this value looks reliable. In any case, these SRC data were used mainly to constrain the opposition surge. We tested the reliability of the derived opposition surge parameters by changing the normalization value to estimate the radiance factor, finding that this factor does not affect the value of B₀ and h parameters (changes are within the error bars).

The amplitude and width of the opposition effect are 2.283 ± 0.029 and 0.05728 ± 0.0005, respectively. The amplitude of the opposition surge is considerably lower than the one reported by Simonelli et al. (1998) from the analysis of global-averaged photometry from the Viking mission data acquired with the clear filter in the 1.5–123° phase angle range ($B_0=4_{-1}^{+6}$, h = 0.055 ± 0.025). These results were determined with very few observations at small phase angles, one observation at α = 1.5° and a couple at α = 3°, thus originating the large error bars. In the Hapke theory, B₀ was originally expected to vary between zero and one (Hapke 1986), and later models relaxed this constraint to the zero-to-three range, especially for bodies characterized by irregular surfaces, as commonly found for several asteroids, comets, and planets satellites (Simonelli et al. 1998; Li et al. 2015; Fornasier et al. 2015).

The SRC data are unique in the coverage of Phobos’ small phase angles, and even if not absolutely calibrated, they provide the best estimation of the amplitude and width of Phobos’ opposition surge so far. The parameters of the opposition surge of Phobos are very similar to the ones found for comet 67P from disk-resolved photometry using the Hapke (2012) model (Fornasier et al. 2015). Cometary nuclei are in fact very red objects, just as Phobos is, and have a comparable albedo, even if their surface is usually even darker (6.5% at 649 nm for comet 67P). For the disk-averaged photometry of 67P comet, Fornasier et al. (2015) found an opposition amplitude ranging from 1.91–2.22 in the 480–989 nm range but a width (0.021–0.032) smaller than the one we report for Phobos.

In Fig. 3, we show the observations at the opposition as well as some SRC data acquired at higher phase angles to test the goodness of the model. These data have been corrected for the exposure time and the square of the heliocentric distance and normalized to match the K076 observations at α = 16–18°. Even if the data are scattered, the Hapke fit appears to be good enough, even at high phase angles.

To model the HRSC data, we fixed the opposition parameters to the values found with the SRC camera, and we computed the SSA and the asymmetric factor. We also made an attempt to model the average roughness slope; however, the Mars Express data cover sixteen years of observations and include very different geometry and varying stray light contributions from Mars, resulting in scattering of the data and very large uncertainties in the roughness parameter. Therefore, we decided to fix θ–24°, the value determined for the green filter from the disk-resolved photometry. This value is also within the error bars of the roughness determined by Simonelli et al. (1998) from the Viking data. The results of the modeling are reported in Table 1 and shown in Fig. 4.

Besides the strong opposition effect, Phobos shows a relatively dark and backscattering surface. The SSA ranges from 6.4% in the blue filter to 9.1% in the IR one, increasing with the wavelength, a behavior expected because of Phobos’ red spectrum. The asymmetric factor 𝑔 is negative, indicating backscattering of the incoming light, and slightly decreases at higher wavelengths. The SSA in the green filter is 7.4%. We compared the values of the reflectance in the green filter with those published by Pajola et al. (2013) from observations acquired with the OSIRIS cameras on board Rosetta. Their green filter data was centered at the same wavelength as that of Mars Express but narrower, and the authors reported a reflectance value of 0.028 at α = ~19°, which is very similar to the one we report here that was derived at a similar phase angle (Fig. 4). It should be noted that the region observed by Rosetta is in the trailing side and is not observable by Mars Express because it orbits between Mars and Phobos (Mars Express observations cover mostly the leading side of the satellite). Therefore, these similar values indicate that, on a global scale, the reflectance properties of the red unit are very close in the trailing and leading sides. When comparing these parameters with those obtained from the disk-averaged photometry of other dark solar system bodies (Table 1), Phobos appears quite peculiar, having a stronger opposition effect than the C-type asteroids investigated, excluding Mathilde, for which very few data were available to constrain the opposition surge, and therefore it has large uncertainties in B₀ (Clark et al. 1999). Its steep opposition surge is similar to that of comet 67P and is likely due to a very porous surface that favors a higher shadowing than in dark asteroids. The SSA is higher than that of cometary nuclei or C-type asteroids but comparable to that of the dwarf planet Ceres. The photometric parameters of the distinct Mars terrains reported in the literature are very different from those of Phobos. Mars terrains have much higher SSA values and are less backscattering, and they usually have a lower roughness and a smaller opposition surge (Johnson et al. 2006a,b; Jehl et al. 2008).

Fig. 3 Integrated phase function of Phobos from SRC observations. Top: opposition effect from integrated photometry from the five SRC orbits. The data presented in Fig. 2 have been normalized, and the radiance factor was estimated from the measurement with the green filters at α = 5°. Bottom: SRC phase function (black points) superposed on the Hapke global fit model. For completeness, some SRC data acquired during 2019–2020 at a higher phase angle are shown on the gray background, but these data were not used in the analysis.

Fig. 4 Phobos phase curve for the HRSC observations (black points) acquired with the blue, green, red, and IR filters. The red line represents the Hapke global fit model.

Fig. 5 Phobos phase curve in the blue (top-left), green (top-right), red (bottom-left), and IR (bottom-right) filters, respectively. The continuous line represents the HG IAU model that best fits the data, while the blue dashed line represents the linear fit to the data for phase angle >7º.

4 Absolute magnitude

To evaluate the absolute magnitude of Phobos in the four calibrated filters, we first computed the magnitude of the sun in the HRSC filter by convolution of the solar spectrum with the filter transmission. Then, we corrected the flux considering the fact that Phobos’ spectrum is redder than that of the Sun using $$F_{\text{corr }}=F_{\text{unc }}\times \frac{{\displaystyle {\int }_{\lambda }{F}_{\odot }}(\lambda )T(\lambda )\text{d}\lambda }{{\displaystyle {\int }_{\lambda }{F}_{\text{phobos }}}(\lambda )T(\lambda )\text{d}\lambda },$$(2)

where F_unc and F_corr are respectively the uncorrected and corrected Phobos reflectance factors at the central wavelength λ_c of the considered filter, T(λ) is the filter throughput, and F_phobos(λ) and F_⊙(λ) are the Phobos red unit (Fraeman et al. 2012) and the solar (from the HST catalog) spectra, respectively, both normalized to unity at the central wavelength of the considered filter (λ_c).

We computed the magnitude of Phobos in the given HRSC filter as $$m_{\text{filt }}=-2.5\text{\hspace{0.17em}log}\frac{\left(F_{\text{corr }}{\displaystyle {\int }_{\lambda }{F}_{\text{phobos }}}(\lambda )T_V(\lambda )\text{d}\lambda \right)}{\left({\displaystyle {\int }_{\lambda }{F}_{\text{vega }}}(\lambda )T_V(\lambda )\text{d}\lambda \right)}\text{, }$$(3)

where F_vega is the flux of Vega from the HST catalog.

The reduced absolute magnitudes in the four filters are shown in Fig. 5. We characterized the phase function behavior using the IAU HG system (Bowell et al. 1989), and the associated results are reported in Table 2. We also report the linear magnitude (H_lin), that is, the magnitude computed neglecting the opposition surge; the linear slope coefficient (β), calculated for α > 7°; and the estimated phase integral q using the equation of Bowell et al. (1989): $$q=0.29+0.684\times G.$$(4)

The G parameter for Phobos ranges from 0.025 for the IR filter to 0.032 for the blue filter, the averaged value of the phase integral estimated from the G value is 0.308, and the slope parameter β is 0.0307 ± 0.001 mag/°. The slope parameter and the phase integral are similar to the values reported in the literature for primordial C-type asteroids. Vereš et al. (2015) found on average a β value of ~0.03 mag/°for C-type asteroids, while Shevchenko et al. (2019) obtained a phase integral ranging from 0.031 to 0.037 for dark asteroids. The phase integral that we determined is very close to the value reported for Phobos by Simonelli et al. (1998) (q = 0.30 ± 0.04), while Shevchenko et al. (2019) found a higher value (q = 0.38 ± 0.03) using the HG1G2 system. The slope parameter is similar to what is already reported in the literature (for instance β = 0.032 and 0.033 mag/deg reported by Noland & Veverka 1976 and Simonelli et al. 1998, respectively). The G value for the primordial asteroids varies depending on the sources. On average, it is 0.09 for D-type asteroids (Vereš et al. 2015), while for C-type asteroids it ranges from 0.07 (Shevchenko & Lupishko 1998) to 0.15 (Vereš et al. 2015), and Ciarniello et al. (2017) reported a G value of 0.02 for Ceres. Phobos thus has a moderate brightness dependence on the phase angle; this dependence is less steep than the one of comets, which have very similar spectral behavior in the visible range. For instance, comet 67P/Churyumov-Gerasimenko has a negative slope value (−0.13 ± 0.01) and a steeper phase coefficient (β = 0.047 ± 0.002 mag/°, Fornasier et al. 2015).

5 Disk-resolved photometry with HRSC

The calibrated disk-resolved data obtained from the HRSC camera were studied using the latest Hapke IMSA model (Hapke 2012). The photometric model was used to characterize the surface through a set of photometric parameters. We created cubes of data containing the original image, phase, emission, incidence, latitude, and longitude images after they had been co-registered following the procedure explained in Sect. 2.1. We filtered the data for each cube, excluding the pixels with an incidence and/or emission angle greater than 70° in order to avoid unfavorable observation geometries near the limb or the terminator. We also excluded a number of images for which the Phobos shape in the corresponding simulations did not reproduce the original one accurately enough as well as images where the co-registration between the original image and the simulated one was not satisfactory.

The Hapke IMSA model (Hapke 2012) allows for the estimation of the porosity of the uppermost layers via the porosity factor K defined as: $$K=-\ln \left(1-1.209{\varphi }^{2/3}\right)/\left(1.209{\varphi }^{2/3}\right),$$(5)

where ϕ is the filling factor and the surface porosity is defined as 1-ϕ (Helfenstein & Shepard 2011). In the model, the opposition effect is reproduced considering both the shadow hiding and the coherent backscattering, each having two parameters describing the width and the amplitude of the opposition surge. However, the calibrated data poorly cover the small phase angle range, and it is not possible to correctly evaluate it. For this reason, we constrained the opposition surge using the parameters derived from the SRC camera disk-averaged photometry, therefore considering only the shadow-hiding effect. Neglecting the coherent backscattering implies underestimating the opposition surge at very small phase angles (<2°). However, it is well known that multiple-scattering is very low for dark surfaces, such as the one of Phobos; therefore, the coherent backscattering should have a negligible contribution in the opposition surge of Phobos. We modeled all calibrated filter images using a least- χ² fit based on the Levenberg-Marquardt algorithm. We looked for solutions within the following boundaries: w_λ = {0.01, 0.3}, 𝑔_λ = {−1.0, 1.0}, and $\overline{\theta }=\left\{5°,90°\right\}$.

The uncertainties of the Hapke parameters in the disk-resolved case were determined using a bootstrap method (Efron & Tibshirani 1993). In fact, the errors determined from the covariance matrix were extremely small because of the large number of points (> 10 million) fitted for each filter. In addition, due to the degeneracy in the Hapke model parameters, the fit sometimes converges to a local minimum instead of a global minimum. Therefore, to avoid these problems, we resampled the data, generating smaller sets each including thousands of data in reflectance that we fitted using the Hapke (2012) model. We also resized the boundaries of the parameters closer to the solution found, that is, w₀ = 0.06 ± 0.04, 𝑔₀ = −0.269 ± 0.1, and θ₀ = 24 ± 8°, and we slightly varied the input parameters for each run of modeling. The resampling was performed ten times, resulting in 100 runs of the fit routine with different samples and initial conditions (examples of the bootstrap data selection and modeling are reported in Fig A.2). The modeling for the different subsets of resampled data gave stable solutions, and finally the parameter values and uncertainties were determined from the mean and the standard deviation of 100 iterations. The best solutions for Phobos disk-resolved photometry are reported in Table 3 and shown in Fig. 6.

Similar to what was found for the disk-averaged solution, the w_λ increases with the wavelength, though in a less steep way. The asymmetric factor indicates backscattering, and its value is almost constant in the different wavelengths, considering the error bars. The roughness parameter varies between 22° and 26.4°, but with relatively large error bars; therefore, it may be considered almost constant with the wavelength, as expected because it is a geometric parameter in the model and therefore not wavelength dependent.

Phobos also has a very high surface porosity (87%), a value very similar to the one reported in the literature for comet 67P (Fornasier et al. 2015), indicating that Phobos’ top surface is covered by a thick dust mantle with grains that likely have a complex structure, such as fractal aggregates.

Fig. 6 Resolved photometry (gray points) of Phobos with the four calibrated HRSC filters, where each point represents the reflectance of an element of the Phobos surface at a given incidence, emission, and phase angle, with the modeled data superposed in red (Hapke 2012). These simulations were performed using the 3D shape model (Ernst et al. 2023) and the best-fit solution of the Hapke modeling for each filter (Table 3).

5.1 Single-scattering albedo maps: Phobos reflectance variations

In the next step, SSA maps were produced in all four HRSC filters. The SSA maps were created by reversing the Hapke (2012) model equation: $$\begin{array}{l}\omega =\frac{I/F}{\frac{K}{4}\cdot S(i,e,\alpha ,\overline{\theta })}\hfill \\ \times \frac{1}{P_{\text{hg}}(\alpha ,\omega )\cdot B_{\text{sh}}\left(\alpha ,B_{\text{sh},0},h_{\text{sh}}\right)+B_{\text{cb}}\left(\alpha ,B_{\text{cb},0},h_{\text{cb}}\right)\cdot M\left(\frac{{\mu }_0}{K},\frac{\mu }{K},\omega \right)}.\hfill \end{array}$$(6)

The set of photometric parameters are from the solutions of the resolved photometric analysis previously determined in each filter. To generate the SSA maps, we carefully selected the images, ensuring good spatial resolution and phase angles that were not too large in order to avoid the important shadowing effect in the analysis. It is worth noting that the trailing side is only partially visible from Mars Express due to the fact that this spacecraft is in a polar orbit positioned between Mars and Phobos. The images used for the maps are listed in Tables B.1–B.4 available at CDS. For overlapping areas observed in more than one image in a given filter, we selected the best resolved image and cropped it accordingly. The blue and green filters have the best image quality in terms of resolution and surface coverage compared to the red and IR filters. The final SSA maps are presented in Figs. 7, and A.3.

The SSA maps correctly display the main geomorphological features of Phobos, including the Stickney and Limtoc craters and the Kepler dorsum (Fig. 7). Considering the map generated with the green filter, the average SSA (ω) is 7.1%, its value ranges from 5.5 to 8% for the red unit, and the blue unit is distinctly brighter, as previously observed in the literature (Fraeman et al. 2014; Simonelli et al. 1998; Pajola et al. 2013). Notably, the northeast rim of the Stickney crater exhibits a higher reflectance than any other region on Phobos, with ω higher than 10%. Some crater rims in the west side of the Stickney crater are also brighter, with ω ~ 0.010 in the green filter. The floor of the Stickney crater belongs to the darkest area, with ω ~ 5.5–6%. The regions located on the western part of Stickney are characterized by an SSA value that is between the typical red unit and the darkest region in the Stickney floor (6–6.5%).

All of these surface albedo heterogeneities are also evidenced in the images corrected for the illumination conditions using the Lommel-Seeliger law and without applying any phase angle corrections (Fig. 8).

To better study local reflectance heterogeneities, we also considered the SRC data in SSA relative to the mean value of 7.1%, which is the average ω in the green filter. Figure 9 maps the relative ω variations from one image collected during the K076 orbit in November 2019 at a spatial resolution of 22 m px⁻¹ and one from orbit L155 acquired in September 2020 at a resolution of 32 m px⁻¹. Comparing these maps with the main geomorphological features and color units reported in Basilevsky et al. (2014), we confirm that the areas dominated by the blue unit have a higher reflectance. The brightest areas of Phobos, with ω 50–60% higher than average, can be found in the northeast Stickney rim, where the blue unit dominates, and observed in relatively small craters producing blue material ejecta as well as on the rims of some grooves. The northwest rim of the Stickney and Limtoc craters also are among the brightest areas, with ω being locally 40–50% brighter than average. The south rim of the Drunlo crater, the rims of the Reldresal crater, and the blue unit near the Kepler dorsum are moderately brighter, with ω 15–30% higher than average, while the regions south of Stickney are 10–20% brighter. Darker areas can be found near the west side of Stickney, in the equatorial region, where ω is 10–20% darker than average, and in the hummocky terrain characterizing the floor of Stickney, which is on average 10% darker and associated with landslides (Basilevsky et al. 2014) but shows locally brighter spots. Grooves usually have a darker floor and brighter rims. More specifically, the grooves located in the blue unit display brighter southern rims.

Fig. 7 Single-scattering albedo map for the green filter.

Fig. 8 Examples of Phobos’ reflectance factor for areas observed at different phase angles. The I/F value has been corrected for the illumination conditions using the Lommel-Seeliger law.

Fig. 9 SRC map showing the relative variation of the albedo compared to the average value (7.1%). Top: relative albedo variation from the image 2019-11-17T03:27:01 acquired during K076 orbit. Bottom: relative albedo variation from the image 2020-09-26T02:18:40 acquired during L155 orbit.

5.2 Photometric properties of blue and red units

We also investigated the variations in the opposition effect among three different regions of interest (ROIs) on Phobos, two representing the blue unit (northeast rim of Stickney and the brightest rim of the Limtoc crater) and one in the red unit located in the floor of Stickney crater (Fig. A.4). We first analyzed the SRC images since they provide the best coverage of the Phobos opposition effect, especially around the Stickney crater, where the images are particularly well-resolved (spatial resolution of 23 m px⁻¹) at low phase. The K076 orbit is especially noteworthy due to its phase angle coverage, ranging from 1 to 17°. Since these data are not absolutely calibrated, we applied the same calibration factor used for SRC disk-integrated analysis. Because we only have data at small phase angles, we set $\overline{\theta }$ to 24° and g to −0.267 in the SRC ROI-resolved photometric analysis. The asymmetry parameter value was chosen based on the results of the disk-resolved Hapke fitting for the green filter. The results of the SRC disk-resolved analysis (Table 4) confirmed the main findings previously discussed: the Stickney crater’s northeast rim covered by blue unit material is significantly brighter (ω = 0.092) compared to the crater’s floor (ω = 0.059), resulting in an ω difference of ~56%. Apart from the reflectance differences, the phase functions of the different ROIs of Phobos (and consequently of the different Phobos units) exhibit diverse behaviors. In particular, the B₀ parameter appears slightly lower for the blue unit (Stickney rim, B₀ = 2.02) than for the red unit (Stickney floor, B₀ = 2.65). Therefore, the opposition effect is more pronounced for the red unit than for the blue one. Simonelli et al. (1998) previously also found global changes in the phase function among the two units.

These authors also showed that the reflectance of the dark material usually located in the floor of craters tends to darken faster with increasing phase angle than the bright material. The dark material found in the Stickney crater floor and in other smaller craters (Goguen et al. 1978), located mostly in the trailing side and not observable by Mars Express, have been reported to have lower ω, to be more backscattered, and to be characterized by a higher roughness than Phobos’ average terrain (Goguen et al. 1978; Simonelli et al. 1998). We investigated the asymmetric factor and roughness of the three ROIs using the green filter observations. We fixed the opposition parameters to the values determined by SRC observations for each individual ROI, and we ran the Hapke (2012) model, this time leaving ω, 𝑔, and $\overline{\theta }$ as free parameters. Indeed, we found that the SSA at 538 nm is lower in the red unit material composing the floor of Stickney. However, the ROIs located in the blue and red units have the same asymmetric factor value, indicating very similar backscattering properties (Table 4). Additionally, we did not observe differences in the surface porosity for the red and blue unit ROIs investigated in this work. The porosity has in fact a high value indistinguishable from the Phobos’ average one reported in Table 3. Concerning the roughness, our results do not support the conclusions of Goguen et al. (1978) and Simonelli et al. (1998) that the floor of craters is significantly rougher than the average terrain, at least for the floor of Stickney. In fact, our analysis shows that it has a roughness parameter of 18.2°, which is slightly lower but comparable within the uncertainties to the value found for the Stickney rim and is significantly lower than the roughness found in the bright landslide of the Limtoc crater (Table 4). Stickney and Limtoc are among the oldest regions on Phobos, with ages ranging from 2.6 Gyr, if Phobos is a captured asteroid, to 4.3 Gyr, if Phobos was formed in the present orbit (Schmedemann et al. 2014). The ROIs representing the blue unit we investigated include the Stickney ejecta and the bright landslide in Limtoc, therefore potentially relatively younger areas. Our analysis indicates that the blue and red units not only have different spectral behaviors and brightness, but they also have different opposition parameters, with the darker material in the floor of Stickney having a steeper opposition effect than the blue unit ROIs but sharing a similar backscattering behavior. The Stickney floor, resulting from accumulation of downslope material movement, appears to be dominated by darker particles and is smoother than the average Phobos terrain. We note that this area also includes brighter spots of blue material not fully resolved in the images used in this analysis.

6 Discussion

Phobos shows important heterogeneity represented not only by the two distinct spectral units (blue and red) but also in the albedo. The brightest regions, generally found on Stickney’s northeast rim, can have a reflectance 50–65% higher than the average terrain, while darker material is mostly found inside craters and in the regions located at midlatitude and on the dark red unit located at the west side of Stickney.

The selected ROIs within the blue and red spectral units display distinct photometric characteristics. Specifically, the red unit exhibits a higher opposition surge and a lower reflectance compared to the brighter blue unit. In what follows, we compare the Stickney crater properties with those of craters on other small bodies observed by space missions. The Phobos surface is under the gravity dominated regime (diameter ≳ 10 km), where cratering generally produces roundish craters with well-defined ejecta blankets and rims (Holsapple 1993). Other objects under the same regime and recently explored by space missions are (2867) Steins (ESA/Rosetta), (21) Lutetia (ESA/Rosetta), (4) Vesta (NASA/DAWN), and (1) Ceres (NASA/DAWN) and, in the Kuiper belt, the classical Kuiper belt object (486958) Arrokoth, observed by NASA’s New Horizons spacecraft. Even though these bodies are very different in composition, formation processes, and evolution, we try to qualitatively compare the photometric properties of the craters observed on them.

From the aforementioned sample, only Ceres and Arrokoth have primitive carbonaceous organic-rich analog compositions. Ceres is a 965-km dwarf planet whose surface is scattered with very defined craters. Occator is the most intriguing crater on the Ceres surface. The bright central faculae connected to brines and fresher material make this crater particularly unique. On the other hand, Narwish and Kerwan quadrangles are considered to be much older (~2 Ga, Carrozzo et al. 2019 and ~1 Ga, Williams et al. 2018). From the maps of the photometric parameters produced by Li et al. (2019), the Narwish region is characterized by a lower w, stronger 𝑔, and higher $\overline{\theta }$ than Ceres. Kerwan also displays a high $\overline{\theta }$ but lower 𝑔 and higher w. Narwish’s photometric parameters are similar to the characteristics of the Stickney floor, where w is proportionally lower and B₀ is higher. While B₀ and $\overline{\theta }$ do not probe the same size scale, a higher B0 value in dark bodies is generally related to rougher and irregular particles and to the increase of mutual inter-particle shadows.

For (486958) Arrokoth, which has approximately a size similar to that of Phobos but a complex bilobated shape, observations are available in only three different phase angles (12°, 16°, and 33°), thus limiting the comparison to only albedo variation. The largest noticeable crater is found in the smallest lobe and is called Sky. Sky shows two apparent landslides at its crater walls, landslides which are 60% brighter than the overall Arrokoth (Hofgartner et al. 2021). The crater floor is partially hidden, but it is apparently neutral or slightly darker. Phobos also shows bright landslides in the eastern walls of Stickney and Limtoc, reaching 40–60% higher reflectance than the average terrain. Landslides in crater walls are apparently ubiquitous, and they are also present in (21) Lutetia, (4) Vesta, and Ceres, but they are darker in the Lutetia Baetica region (Hasselmann et al. 2016) and varied in Ceres as well as Vesta. Dark landslides have been proposed to be an outcome of particle sorting (Hasselmann et al. 2016): small particles are levitated during the event, thus becoming lost or transported away. Bright landslides, on the other hand, are a possible outcome of fresher materials or diverse subsurface composition becoming exposed, as observed on Phobos.

We also attempted to compare Phobos’ photometric properties with those of meteorites and minerals. The bidirectional distribution functions (BRDF) measured in the laboratory on different materials, grain sizes, and geometric conditions demonstrate that it is difficult to obtain a phase function similar to that observed on Phobos. In fact, the scattering function is usually more isotropic in the laboratory (Souchon et al. 2011; Johnson et al. 2013; Potin et al. 2022; Wargnier et al. 2023), resulting in a significant increase in reflectance after a 80–100° phase angle due to forward scattering. Moreover, laboratory phase function measurements indicate that the opposition surge is often less significant than expected or observed for Solar System bodies. For instance, Beck et al. (2012) found a B₀ parameter of approximately 0.4 for the Tagish Lake meteorite, whereas observations of primitive asteroids, from which the meteorite could have originated, typically exhibit values greater than 1.5. None of the different meteorites measured by Beck et al. (2012) (including carbonaceous chondrites, Tagish Lake meteorite, howardites, eucrites, diogenites, and a lunar sample) show an opposition surge as strong as the one observed for Phobos. If we computed the absolute and relative opposition effect intensity (OEI) from the SRC observations of Phobos using the same definition as Beck et al. (2012)³, we obtain for Phobos a value of ~0.06 for the absolute OEI and ~3.7 for the relative OEI, values that are much higher than what was reported for the meteorites measured by Beck et al. (2012) (see their Fig. 6). It should be noted that the reflectance at zero phase angle is often extrapolated in laboratory data, where measurements usually cannot be performed at phase angles smaller than ~5° and different parameters such as grain size and shape, surface rugosity, and porosity affect the shape of the phase function (Okada et al. 2006; Kamei & Nakamura 2002). For instance, studies conducted on different grain sizes of the same sample have shown that the forward scattering tends to be larger for the finer grains (Kamei & Nakamura 2002; Johnson et al. 2013), and that larger particles tend to diffuse more backward (Souchon et al. 2011). Moreover, for smaller grains, the effect of surface charging on particle scattering can also induce discrepancies between laboratory and remote sensing data (Klačka & Kocifaj 2007).

Overall, our photometric analysis has shown that Phobos photometric properties show a close resemblance to those of comet 67P: both have a red spectrum, a high surface porosity, and similar opposition effect values, even though it should be noted that the material on top of comet 67P has a lower SSA (3.4–4.2%) compared to Phobos (7.4%), and it is more backscattering (𝑔 = −0.42 for 67P, and −0.27 for Phobos). As discussed in the previous sections, the SSA of Phobos is higher than other primordial bodies investigated by space missions (Bennu, Ryugu, Mathilde) and comparable to that of the dwarf planet Ceres.

Figure 10 shows Phobos’ reflectance ratio at a small phase angle versus the geometric albedo compared to the properties of different asteroidal classes (Belskaya & Shevchenko 2000) and comet 67P (Ciarniello et al. 2015; Fornasier et al. 2015). Phobos has a strong opposition surge that is usually higher than most carbonaceous asteroids investigated and is comparable in intensity to that of metallic or silicate-rich asteroids. However, its albedo is much lower than the one of silicate or metallic asteroids and similar to that of primordial asteroids. Concerning D-type asteroids, which are characterized by a red spectrum (as also observed for Phobos), there are very few data in the literature. Both Shevchenko et al. (2012) and Mottola et al. (2023) reported observations at a small phase angle for Jupiter Trojans (including four D-types), finding that these bodies have a linear dependence in the phase function with almost no opposition surge and usually a very low albedo.

Actually, comet 67P is the body closest to Phobos in terms of photometric properties when considering both the albedo and opposition surge properties. This similarity is also confirmed when comparing the geometric albedo and the phase integral of different dark Solar System bodies (Fig. 10). In this paper, we have derived a phase integral of 0.23 in the green filter from Hapke-resolved photometry and of 0.31 from the HG fitting of magnitude in the same filter from integrated photometry. In the literature, there are several phase integral estimations of Phobos with quite different values: Pang et al. (1983) reported q = 0.33 and Klaasen & Duxbury (1979) q = 0.27, while Shevchenko et al. (2019) reported q = 0.38 ± 0.03 (from HG1G2 model) and Simonelli et al. (1998) q = 0.3. Deimos has a phase integral of 0.38 (Thomas et al. 1996).

Among the two values reported in our analysis, we think that the one derived from the Hapke-resolved photometry should be more reliable because it is based on data acquired at a very small phase angle, therefore permitting accurate modeling of the Phobos strong opposition effect, while the HG model was based on the reduced magnitudes derived from absolutely calibrated filters, which have a limited coverage of the opposition surge. When comparing the phase integral and geometric albedo of Phobos with those of other dark Solar System bodies (Fig. 10), it appears once more that Phobos has photometric properties similar to those of Jupiter family comets, which originate from the Kuiper belt. Phobos’ phase integral is also very close to that of Phoebe, which is supposed to be a Kuiper belt object captured by Saturn. Additionally, we also found that the opposition effect of Phobos is very similar in shape and parameters to the one of comet 67P, and it has a similar high surface porosity. In Fig. 11 we also compare Phobos’ red unit spectral slope, derived from Fraeman et al. (2012) spectra, to that of other primordial bodies, derived from the analysis of Fornasier et al. (2007), such as Jupiter Trojans, which are dominated by D-type asteroids; cometary nuclei; and the different dynamical classes of the transneptunian population. The blue star refers to the spectral slope of comet 67P determined at ~3.6 au pre-perihelion and at a phase angle of 1.3° (Fornasier et al. 2015). As already well known, Phobos has a red spectral behavior close to that of D-type asteroids, and it also shares spectral similarities with cometary nuclei, notably comet 67P; with Jupiter Trojans; and with the moderately red bodies of the transneptunian population.

Based on the aforementioned similarities in the photometric and spectral properties with primordial Solar System bodies, particularly with cometary nuclei, we suggest that additional dynamical investigations deserve to be performed in order to understand the origin of Mars satellites, mimicking the capture not only of asteroids but of other small bodies. These simulations may consider a bilobated or binary comet or external Solar System body capture to test the hypothesis that Phobos and Deimos may currently be exhausted cometary nuclei captured by Mars. A binary or bilobated comet might eventually be fragmented into two bodies by Mars’ tidal and gravitational forces during the capture process and may have subsequently generated the Martian moons.

Fig. 10 Comparison of Phobos photometric properties with those of other Solar System bodies. Top: Phobos reflectance ratio (at α = 0.3º over α = 5º) versus geometric albedo compared to the values observed for different classes of asteroids (following the Tholen & Barucci 1989 classification scheme), derived from Belskaya & Shevchenko (2000), and comet 67P (Fornasier et al. 2015). Bottom: figure adapted from Verbiscer et al. (2022) showing the phase integral versus geometric albedo for dark Solar System bodies. The Phobos value derived in this paper from the green filter and using the Hapke modeling of resolved photometry is represented by a red star (uncertainties are within the star symbol). Literature data of Phobos (q = 0.30) from Simonelli et al. (1998) are also reported. We refer to Verbiscer et al. (2022) (their Fig. 11) for the references of the individual phase integral and geometric albedo values of the objects.

Fig. 11 Spectral gradient, estimated in the 550–800 nm range, of Phobos’ red unit (11%/(100 nm); red vertical line) derived from spectra presented in Fraeman et al. (2012) compared to the values reported by Fornasier et al. (2007) for cometary nuclei, Jupiter Trojans, and the different dynamical classes of the transneptunian population. Frequency refers to the number of objects in the different histograms. The figure was adapted from Fornasier et al. (2007). The blue star indicates the spectral slope of comet 67P determined at phase 1.3º (12.6%/(100 nm)) and at a heliocentric distance of 3.6 au by Fornasier et al. (2015).

7 Conclusions

We have analyzed photometric observations performed by the Mars Express mission on Phobos with the HRSC absolutely calibrated filters in the blue, green, red, and IR range and with relatively calibrated panchromatic data acquired with the SRC camera. In this work, we have presented the results of the photometric modeling for the different HRSC filters, the SSA maps, and a comparison of the photometric properties of Phobos with those of other satellites and dark minor bodies. The data acquired with the SRC camera are unique and represent the best coverage achieved so far of the Phobos opposition effect. This camera in fact acquired 299 images covering the 0.3–17° phase range. The main findings of our work are the following:

• Phobos has a strong opposition effect characterized by an amplitude B₀ = 2.283 ± 0.029 and a width 0.05728 ± 0.0005.

• Phobos is a relatively dark object with a mean albedo in the green filter of 6.83%. The SSA increases with wavelength from 6.2% in the blue filter to 8% in the IR filter, as expected because of Phobos’ very red spectrum.

• Its surface is backscattering (−0.027 < 𝑔 < = −0.024) and has a high porosity value (87%), indicating the presence of a thick dust mantle, possibly composed of grains with a complex shape or fractal aggregates.

• The SSA maps reproduced the main geomorphological features of the Mars satellite and provide clear evidence of the albedo dichotomy between the blue unit, which is up to 50–65% brighter than average in the northeast rim of Stickney, and the darker red unit. The darkest regions are found in the floor of Stickney and in the regions located in the western side of this crater.

• Local photometric analysis of selected areas located in the blue and red units showed that the red unit terrain has a lower SSA and a stronger opposition surge compared to the blue unit terrain, but they have similar surface porosity values and backscattering properties.

• The phase function of Phobos shares analogies with dark asteroids, but its SSA is considerably higher than the values reported for carbonaceous-rich asteroids, such as Bennu, Ryugu, and Mathilde, and similar to the value found for the dwarf planet Ceres.

• The opposition surge parameters, the porosity, and the phase integral of Phobos are very close to the value reported in the literature for the comet 67P, which is characterized by a similarly red spectral slope but has a much lower SSA value.

• The phase integral of Phobos is similar to that of dark asteroids of the Jupiter family comets and Phoebe. While waiting for detailed in situ analysis of D-type Jupiter Trojans by the Lucy mission, so far the closest analogs to Phobos photometric properties are cometary nuclei. However, the material composing Phobos has a higher SSA than that on the top surface of comets or primordial asteroids.

• Based on the photometric properties that are similar to the comet 67P, we suggest investigating the dynamical capture of a binary or bilobated comet as the potential origin of the Martian moons.

All of these results are of high interest and support to the JAXA MMX mission. The mission’s main goals are the return of Phobos samples collected in both the red and blue units, a detailed investigation of Mars’ satellites, and determining the origin of the Martian moons.

Acknowledgements

We acknowledge the ESA Planetary Science Archive for space mission data procurement and the principal investigator of the HRSC instrument, G. Neukum (Freie Universitaet, Berlin, Germany). We thank the anonymous reviewer for valuable comments and suggestions, which helped us to improve the quality of the manuscript.

Appendix A Supplementary material: Figures

Fig. A.1 The image on the left/right is an example of original/simulated image (incidence angle). The simulated image is not satisfactory. In fact, Phobos’ shape is larger in the simulated image than in the original one, and it was thus discarded from our analysis.

Fig. A.2 Example of bootstrap method for the green filter resolved photometry. The whole sample, which included several million points, was resampled to generate smaller datasets of approximately one thousand points each. Here, two subsets are shown (gray points) together with the best-fit data models (red points).

Fig. A.3 Single-scattering albedo maps in the blue (top), red (center), and IR (bottom) filters.

Fig. A.4 Location of the three ROIs for which the local resolved photometric properties have been investigated (see Table 4).

Fig. A.5 Maps generated from SRC images: Relative variation of the albedo compared to the average value (7.1%) determined with the green filter and used to normalize the SRC data. Top: From the image 2019-11-17T03:27:01 acquired during K076 orbit. Bottom: From the image 2020-09-26T02:18:40 acquired during L155 orbit. The gray color scale permits better visualization of the albedo-geomorphological feature correlation compared to the colored map presented in Fig. 9.

References

All Tables

All Figures

	Fig. 1 Images showing the flux increasing with a decreasing phase angle from images acquired during the orbit K076 (top panel) and M 247 (bottom panel) with SRC. The flux is in DU s−1; no other normalization was applied here (the heliocentric distance was the same for a given orbit).
In the text

	Fig. 2 Observations of Phobos opposition effect from five individual SRC observations. Data have been normalized by the exposure time and corrected for heliocentric distance to the square. The insert shows very small phase angle data acquired for orbits L132 and M 247. Data acquired at a phase angle lower than the angular size of the Sun seen from Phobos distance were discarded in the modeling.
In the text

Fig. 3 Integrated phase function of Phobos from SRC observations. Top: opposition effect from integrated photometry from the five SRC orbits. The data presented in Fig. 2 have been normalized, and the radiance factor was estimated from the measurement with the green filters at α = 5°. Bottom: SRC phase function (black points) superposed on the Hapke global fit model. For completeness, some SRC data acquired during 2019–2020 at a higher phase angle are shown on the gray background, but these data were not used in the analysis.

In the text

	Fig. 4 Phobos phase curve for the HRSC observations (black points) acquired with the blue, green, red, and IR filters. The red line represents the Hapke global fit model.
In the text

	Fig. 5 Phobos phase curve in the blue (top-left), green (top-right), red (bottom-left), and IR (bottom-right) filters, respectively. The continuous line represents the HG IAU model that best fits the data, while the blue dashed line represents the linear fit to the data for phase angle >7º.
In the text

Fig. 6 Resolved photometry (gray points) of Phobos with the four calibrated HRSC filters, where each point represents the reflectance of an element of the Phobos surface at a given incidence, emission, and phase angle, with the modeled data superposed in red (Hapke 2012). These simulations were performed using the 3D shape model (Ernst et al. 2023) and the best-fit solution of the Hapke modeling for each filter (Table 3).

In the text

	Fig. 7 Single-scattering albedo map for the green filter.
In the text

	Fig. 8 Examples of Phobos’ reflectance factor for areas observed at different phase angles. The I/F value has been corrected for the illumination conditions using the Lommel-Seeliger law.
In the text

	Fig. 9 SRC map showing the relative variation of the albedo compared to the average value (7.1%). Top: relative albedo variation from the image 2019-11-17T03:27:01 acquired during K076 orbit. Bottom: relative albedo variation from the image 2020-09-26T02:18:40 acquired during L155 orbit.
In the text

Fig. 10 Comparison of Phobos photometric properties with those of other Solar System bodies. Top: Phobos reflectance ratio (at α = 0.3º over α = 5º) versus geometric albedo compared to the values observed for different classes of asteroids (following the Tholen & Barucci 1989 classification scheme), derived from Belskaya & Shevchenko (2000), and comet 67P (Fornasier et al. 2015). Bottom: figure adapted from Verbiscer et al. (2022) showing the phase integral versus geometric albedo for dark Solar System bodies. The Phobos value derived in this paper from the green filter and using the Hapke modeling of resolved photometry is represented by a red star (uncertainties are within the star symbol). Literature data of Phobos (q = 0.30) from Simonelli et al. (1998) are also reported. We refer to Verbiscer et al. (2022) (their Fig. 11) for the references of the individual phase integral and geometric albedo values of the objects.

In the text

Fig. 11 Spectral gradient, estimated in the 550–800 nm range, of Phobos’ red unit (11%/(100 nm); red vertical line) derived from spectra presented in Fraeman et al. (2012) compared to the values reported by Fornasier et al. (2007) for cometary nuclei, Jupiter Trojans, and the different dynamical classes of the transneptunian population. Frequency refers to the number of objects in the different histograms. The figure was adapted from Fornasier et al. (2007). The blue star indicates the spectral slope of comet 67P determined at phase 1.3º (12.6%/(100 nm)) and at a heliocentric distance of 3.6 au by Fornasier et al. (2015).

In the text

	Fig. A.1 The image on the left/right is an example of original/simulated image (incidence angle). The simulated image is not satisfactory. In fact, Phobos’ shape is larger in the simulated image than in the original one, and it was thus discarded from our analysis.
In the text

	Fig. A.2 Example of bootstrap method for the green filter resolved photometry. The whole sample, which included several million points, was resampled to generate smaller datasets of approximately one thousand points each. Here, two subsets are shown (gray points) together with the best-fit data models (red points).
In the text

	Fig. A.3 Single-scattering albedo maps in the blue (top), red (center), and IR (bottom) filters.
In the text

	Fig. A.4 Location of the three ROIs for which the local resolved photometric properties have been investigated (see Table 4).
In the text

Fig. A.5 Maps generated from SRC images: Relative variation of the albedo compared to the average value (7.1%) determined with the green filter and used to normalize the SRC data. Top: From the image 2019-11-17T03:27:01 acquired during K076 orbit. Bottom: From the image 2020-09-26T02:18:40 acquired during L155 orbit. The gray color scale permits better visualization of the albedo-geomorphological feature correlation compared to the colored map presented in Fig. 9.

In the text