MMX geodesy investigations: science requirements and observation strategy

Abstract

In order to investigate the origin of Phobos and Deimos, the Japanese Martian Moons eXploration (MMX) mission is scheduled for launch in 2024. MMX will make comprehensive remote-sensing measurements of both moons and return regolith samples from Phobos to Earth. Geodetic measurements of gravity, shape, and rotation parameter of a body provides constraints on its internal structure reflecting its origin and evolution. Moments of inertia are important parameters to constrain the internal mass distribution, but they have not been well determined for the Martian moons yet. We discuss the mission requirements related to the moments of inertia to detect a potential heterogeneity of the mass distribution inside Phobos. We introduce mission instruments and operational strategies to meet the mission requirements. We present a preliminary imaging strategy from a quasi-satellite orbit for a base shape model that is expected to be created at the early stage of the mission. Geodetic products including ephemeris, gravity field, rotation parameter of Phobos, and spacecraft orbit are of importance not only for the geodetic study, but also for interpreting data from various mission instruments and selecting possible landing sites.

Graphical Abstract

Introduction

The origin of Phobos is still an open issue. It may either be a captured asteroid or have formed from a disk of impact ejecta produced by a giant impact (Craddock 2011). The Martian Moons eXploration (MMX) led by Japan Aerospace Exploration Agency (JAXA) will provide a unique opportunity to answer this question from various aspects including remote-sensing observations and returned samples (Kawakatsu et al. 2019; Kuramoto et al. 2021). MMX is scheduled for launch in 2024 and will make progress in our understanding of planetary system formation as well as primordial material transport between the inner- and the outer-part of the early solar system. The internal structure can possibly provide clues on the formation regime and thus on the origin of Phobos. Although it is not straightforward to determine the origin from the internal structure alone, the internal structure provides hints to the origin and the sequence of evolution. Because of its small size and mass, Phobos could have a very puzzling interior structure (Le Maistre et al. 2019). It is expected that the low density ($\sim$1847 kg/m³, "Bulk density and moments of inertia" section) of Phobos is due to the presence of significant macro-porosity and/or water ice in its interior, in both cases the distribution depends on the formation scenario. The presence or absence of such deviations from a homogenous mass distribution will help us trace back the origin of the satellite.

Geodetic observations of shape, gravity and rotation parameter contribute to detecting the signatures of a heterogeneous interior and providing constraints on internal structure models. We can expect larger signatures from the geodetic observations for larger heterogeneous mass distribution. Some of the key parameters related to the internal structure are the moments of inertia (MOI) which constrain the radial distribution of mass within the interior. MOI are not directly measurable, but can be determined by combining measurements of rotation parameters with measurements of degree-2 gravity coefficients. A rotation parameter that is particularly important for description of Phobos’ rotation is the amplitude of short period forced libration in longitude (e.g., Burmeister et al. 2018). Currently, the MOI are not estimated with sufficient accuracy mainly due to large errors in the gravity coefficients (e.g., Pätzold et al. 2014; Jacobson and Lainey 2014; Yang et al. 2019).

In this paper, we review recent geodetic observations of Phobos, describe geodetic mission requirements related to the MOI using a two-layered Phobos model, and show preliminary observation strategies for the geodetic observations with special emphasis on shape modeling. This paper focuses on observations of Phobos. Observations of Deimos will be discussed possibly later after the design of a flyby trajectory has progressed.

Present state of geodetic observation of Phobos

Shape model

A global shape/topography model can be used to calculate the MOI based on various internal structure models. It may also provide insight into surface processes and regolith evolution (e.g., Miyamoto et al. 2007; Ballouz et al. 2019) and together with gravity models constrain the orbital evolution (Hu et al. 2020). Furthermore, an accurate shape/topography models is also necessary for the photometric correction of images, the identification of potential landing sites having minimal surface hazard, the provision of landmarks for image-based spacecraft orbit determination (Konopliv et al. 2018), and the provision of a reference shape to which an altimeter-measured topographic profile is geometrically fitted to improve spacecraft position (Matsumoto et al. 2020), etc.

Because the shape of a body is essential as described above, much shape and topographic modeling effort on Phobos has been made since explorations of the Martian system started. Willner et al. (2010) introduced a shape model in the form of a digital terrain model (DTM) and a ground control point (GCP) network that was derived from images of the super resolution channel (SRC) of the high-resolution stereo camera (HRSC) on Mars Express (MEX). The GCPs of a planetary body are points of which coordinates are known with respect to the body-fixed rotating frame. The mean accuracy of the GCPs are 21.24 m, 16.19 m, and 16.70 m for XYZ axes of the Phobos body-fixed coordinate, respectively. Willner et al. (2014) then introduced a newer shape model with 100 m/pixel resolution and a set of spherical harmonic coefficients to degree and order 45 estimated from images of the SRC and the framing camera on the Viking Orbiter (VO) with the stereophotogrammetry (SPG) method. Although they did not mention the quantitative accuracy of the 2014 version of the shape model, they stated that the model shows a better representation of Phobos' surface features than their previous model of Willner et al. (2010).

Gaskell (2011) also released another shape model of Phobos through NASA Planetary Data System (PDS). This model was originally reconstructed from 224 VO images and 8 Phobos 2 images in 2006 by the stereophotoclinometry (SPC) method (Gaskell et al. 2008). The model consists of over 3 million triangle plates with the average length of plates of 35.5 m. The data producer gives the RMS residual error of landmarks forming the shape model of 25 m.

Ernst et al. (submitted to Earth Planets and Space) created high-resolution SPC shape models of Phobos and Deimos with images taken by six spacecraft from five missions of Viking Orbiter, Phobos 2, Mars Global Surveyor, Mars Express, and Mars Reconnaissance Orbiter. These models are publicly available in the Small Body Mapping Tool (SBMT, Ernst et al. 2018) distributed at the SBMT website (https://sbmt.jhuapl.edu) and will be archived with the NASA Planetary Data System.

Gravity field

The MEX spacecraft orbits Mars in a highly elliptical and quasi-polar orbit, which provides opportunities to perform close flybys of Phobos every 150 days on average (Witasse et al. 2014). Andert et al. (2010) estimated GM (the product of the Newtonian constant of gravitation G and the body’s mass M) of Phobos from MEX tracking data during its close flyby in 2008 at 275 km distance. Their solution is GM = (0.7127 $\pm$ 0.0021) × 10⁻³ km³/s² with a quoted 1$\sigma$ accuracy of 0.3%. Pätzold et al. (2014) attempted to estimate degree-2 gravity coefficients $C_{20}$ and $C_{22}$ along with GM of Phobos from data during another closer flyby in 2010 at 77-km distance. The updated GM value is (0.7072 $\pm$ 0.0013) × 10⁻³ km³/s² (3$\sigma$) with five times better accuracy than that of the 2008 flyby. However, degree-2 gravity coefficients could not be solved for at sufficient accuracy: $C_{20}$ = − 0.1145 $\pm$ 0.1456 (3$\sigma$) and $C_{22}$ = 0.0025 $\pm$ 0.0056 (3$\sigma$). They explained that the high error of the coefficients is due to strong correlation of the three parameters of GM, $C_{20}$ and $C_{22}$. Yang et al. (2019) analyzed MEX tracking data of two combined flybys in 2010 and 2013 and estimated GM, $C_{20}$, and $C_{22}$. The closest approach distance in 2013 flyby was 59 km. They used the result of Shi et al. (2012), who calculated gravity coefficients based on a shape model and constant density assumption, as a priori information. Their results are GM = (7.0765 $\pm$ 0.0075) × 10⁵ m³ s⁻² (3$\sigma$), $C_{20}$ = − 0.1378 $\pm$ 0.0348 (3$\sigma$) and $C_{22}$ = 0.0166 $\pm$ 0.0153 (3$\sigma$). The corresponding 1$\sigma$ error is 0.04%, 8% and 31%, respectively. The better accuracy for $C_{20}$ seems to be due to the decorrelation of the adjusted parameters owing to the two independent flybys at different geometries and due to the closer flyby distance in 2013 than that in 2010.

The orbit of Phobos is sensitive to its own aspherical gravity field, allowing for the estimation of Phobos’ gravity coefficients from observation of its orbital dynamics. Jacobson and Lainey (2014) applied direct numerical integration method and fitted $C_{20}$ and $C_{22}$ of Phobos together with the Phobos and Deimos ephemeris and GMs using Earth-based astrometric observations since 1877 and spacecraft-based imaging and Doppler tracking data. They emphasize that they did not estimate the amplitude of forced libration of Phobos, but adopted the value derived by Willner et al. (2010). The estimated gravity field coefficients are $C_{20}$ = − 0.1057 $\pm$ 0.0067 and $C_{22}$ = 0.01487 $\pm$ 0.0004. The reported 1$\sigma$ accuracies are about 6% and 3% for $C_{20}$ and $C_{22}$, respectively, but they are estimated with a correlation of 0.96 indicating that the estimates are not independent. Jacobson (2010) pointed out that only one of the three parameters $C_{20}$, $C_{22}$, and libration amplitude can be independently estimated from the observed secular rate of periapsis longitude.

Bulk density and moments of inertia

The mass of Phobos calculated from the GM value of Yang et al. (2019) and $G = \left( {6.67408 \pm 0.00047} \right) \times 10^{ - 11} {\text{ m}}^{{3}} {\text{kg}}^{{ - 1}} {\text{s}}^{{ - 2}}$ (Mohr et al. 2016) is $M = \left( {1.0603 \pm 0.0004} \right) \times 10^{16}\, {\text{kg}}$. The 1$\sigma$ uncertainty of the Phobos’ mass is about 0.04% for which the uncertainty of GM is more dominant than that of G. The bulk density of Phobos is $\rho = \left( {1847 \pm 12} \right)\, {\text{kg/m}}^{{3}}$ by applying the volume $V = \left( {5742 \pm 35} \right){\text{ km}}^{{3}}$ calculated from the shape model of Willner et al. (2014). Currently, the error of the bulk density of about 0.6% is dominated by the error of the volume.

Moments of inertia are of interest in developing models of the satellite's interior structure. They are derived from the amplitude of forced libration in longitude and degree-2 gravity coefficients $\left( {C_{20} , C_{22} } \right)$ of Phobos. If the principal moments of inertia are denoted $A < B < C$, the libration amplitude is expressed as (e.g., Duxbury and Callahan 1989; Willner et al. 2010)

$$\begin{array}{*{20}c} {\theta = \frac{2e}{{\left( {1 - \frac{1}{3\gamma }} \right)}},} \\ \end{array}$$

(1)

where $e = 0.01511$ (Jacobson 2010) is the orbital eccentricity of Phobos and $\gamma$ is the dynamical flattening

$$\begin{array}{*{20}c} {\gamma = \frac{B - A}{C}.} \\ \end{array}$$

(2)

It is not possible to determine the absolute value of $A$, $B$, and $C$ from $\theta$ alone. On the other hand, degree-2 gravity coefficients are related to the MOI, mass $M$, and mean radius $R_{0}$ as follows (e.g., Bills and Rubincam 1995):

$$\begin{array}{*{20}c} {C_{20} = \frac{{\left( {B + A} \right) - 2C}}{{2MR_{0}^{2} }},} \\ \end{array}$$

(3)

$$\begin{array}{*{20}c} {C_{22} = \frac{{\left( {B - A} \right)}}{{4MR_{0}^{2} }}.} \\ \end{array}$$

(4)

Then, by combining $\gamma$ as a function of $\theta$ and ($C_{20} ,C_{22}$), each of dimensionless MOI normalized by $MR_{0}^{2}$ can be derived as follows (e.g., Borderies and Yoder 1990):

$$\begin{array}{*{20}c} {\frac{A}{{MR_{0}^{2} }} = \frac{{\gamma C_{20} + 2\left( {2 - \gamma } \right)C_{22} }}{\gamma },} \\ \end{array}$$

(5)

$$\begin{array}{*{20}c} {\frac{B}{{MR_{0}^{2} }} = \frac{{\gamma C_{20} + 2\left( {2 + \gamma } \right)C_{22} }}{\gamma },} \\ \end{array}$$

(6)

$$\begin{array}{*{20}c} {\frac{C}{{MR_{0}^{2} }} = \frac{{4C_{22} }}{\gamma },} \\ \end{array}$$

(7)

where γ is given by Eq. (1). Because of the large uncertainties of degree-2 gravity coefficients, moments of inertia have not been obtained yet.

Amplitude of forced libration in longitude

The rotation parameters of a planetary body describe the time-varying orientation of the body-fixed rotating reference frame relative to an inertial reference frame. One of the parameters that is of interest with regard to the internal mass distribution of Phobos is the amplitude of forced libration in longitude. Phobos is in synchronous spin–orbit resonance and facing the same side to Mars on average. Superimposed on the uniform mean rotation are oscillations called physical librations that are caused by time-varying gravitational torque that Mars exerts on the dynamical figure of Phobos. The longitudinal libration, an east–west oscillation of the body’s surface about the spin axis, has been observed by spacecraft.

The libration amplitude has been estimated from GCP network analysis using multiple image data. Willner et al. (2010) analyzed image data obtained by the SRC of the HRSC on board MEX and VO framing camera to compute coordinates of 665 GCPs together with the forced libration amplitude of Phobos. They tested various plausible values of the libration amplitude and determined a best fitting value of θ = − 1.2° $\pm$ 0.15° that gives the smallest mean residual of the GCPs. Oberst et al. (2014) updated the results with 813 GCPs and derived a libration amplitude of θ = − 1.09° $\pm$ 0.10°. Burmeister et al. (2018) developed a mathematical framework to directly solve for rotation parameters in the inertial reference frame that is more efficient and faster than the traditional method of scanning the parameter spectrum. They computed a forced libration amplitude of θ = − 1.14° $\pm$ 0.03° together with a GCP network of 685 points.

Lainey et al. (2021) developed new ephemerides of Phobos and Deimos that are fitted to data including 14 years of MEX data. During the fitting process, they solved for the libration amplitude as θ = − 1.09° $\pm$ 0.01°. Note that this dynamical approach adopts Phobos’ gravity field based on the shape model of Willner et al. (2014) and assumes a homogeneous density. The value of Lainey et al. (2021) is close to that of Burmeister et al. (2018), but both values do not agree within their respective error margins.

Mission requirements for MMX geodetic study

Although the libration amplitude of Burmeister et al. (2018) is consistent with that calculated using a Phobos shape model assuming a uniform density or homogenous internal structure (Willner et al. 2014), possible heterogeneous interiors include layered structures with homothetic boundary shape to the surface because such structures give the same dynamical flattening $\gamma$ and hence do not change the libration amplitude. In order to determine the required accuracy of libration amplitude $\theta$ and degree-2 gravity coefficients $\left( {C_{20} , C_{22} } \right)$, we explore a simple two-layer structure that is complementary with models proposed by Le Maistre et al. (2019) (see below). We assume that water ice or porosity is confined to either layer, and calculate how much the MOI deviate from the uniform density case.

Phobos’ bulk density of $\left( {1847 \pm 12} \right)$ kg/m³ (section "Bulk density and moments of inertia") is lower than most of the samples of meteoritical carbonaceous material except for the Tagish Lake meteorite, which requires porosity and/or light elements like water ice. If the low bulk density was explained by water ice, its mass fraction is expected to be 10–35% depending on the rocky material’s grain density (Rosenblatt 2011). If the mass distribution inside Phobos is inhomogeneous, e.g., if water ice is concentrated near the surface or the center, we would observe a deviation of the MOI from the value for a homogenous interior. Here, the MOI differences (∆MOI) with respect to a homogenous case are calculated assuming a simplified tri-axial ellipsoidal figure of Phobos a = 13.03 km, b = 11.40 km, c = 9.14 km (Willner et al. 2014). Because the percentage of the MOI differences are the same for x, y, z directions under the assumption of similar ellipsoidal boundary shape as the surface, we treat ∆MOI as singular value. We assume that water ice is confined either to the outer or to the inner layer and that the rock materials’ density is the same for both layers. The water ice mass fraction is varied between 1 and 30%. The libration amplitude for these models is computed as − 1.14739° and thus consistent with the result of Burmeister et al. (2018). The model-based libration amplitude is computed by discretizing Phobos with 10-m-size cubes and its uncertainty is 10^–5 degrees that is evaluated as relative error with respect to the value computed with 1 m cubes.

Figure 1 shows the MOI differences for the case that the outer layer is composed of rock and water ice. When the outer layer thickness is 10% of the semi-principal axes (labeled as Rb = 0.9R), no more than 14 wt.% of water ice can be contained in the layer and the maximum |∆MOI| is about 9%. When the layer boundary is deeper, more water ice can be contained, but the maximum |∆MOI| is still about 16%. Figure 2 shows the result for the case in which the water ice is confined to the inner layer. If the layer boundary is located below 0.6R, no more than 11 wt.% of water ice can be contained in the inner layer and the maximum |∆MOI| is about 8%. When the layer boundary is shallower, more water ice can be contained, but |∆MOI| peaks at about 17%.

Fig. 1

Moments of inertia differences between homogeneous and inhomogeneous Phobos. The inhomogeneity is assumed to be produced by water ice confined to the outer layer. Different line colors indicate different depths of the layer boundary (red: shallow, blue: deep). Different colors in the circle symbols (color bar below the figure) indicate the density of the outer layer

Fig. 2

Same as Fig. 1, but for the case that water ice is confined to the inner layer. Different line colors indicate different depths of the layer boundary (red: deep, blue: shallow)

We also tested cases in which the porosity is responsible for the low bulk density. We calculated ∆MOI due to inhomogeneous distribution of the porosity using a similar two-layer structure. The results depend on the boundary depth and rock density. We assume that the inner layer is porous which is consistent with a formation scenario of Britt et al. (2002) in which the largest pieces reaccreted first after disruption leaving large void inside. If we limit the computed maximum of the porosity as 60%, the maximum |∆MOI| is about 16% given a boundary position of 0.8R and a rock density of 2700 kg/m³ (Fig. 3).

Fig. 3

Same as Fig. 2, but for the case that the inhomogeneity is produced by porosity in the inner layer. Different line colors indicate different depths of the layer boundary (red: shallow, blue: deep). Different colors in the circle symbols (color bar below the figure) indicate the different porosities of the inner layer (60% maximum)

We set the required accuracy of MOI so that we can detect an inhomogeneous density structure due to water ice of 10% mass fraction that is localized within a relatively small volume of Phobos, i.e., R_b = 0.9R for Fig. 1 and R_b = 0.6R for Fig. 2. Considering |∆MOI| for these situations are 6–8%, the required accuracy for MOI is 3%. Because the MOI are a function of libration amplitude and gravity coefficients (Eqs. 5–7), errors of MOI are described by propagating the errors of libration amplitude ($\sigma \theta$) and gravity coefficients $\left( {\sigma C_{20} , \sigma C_{22} } \right)$ as follows:

$$\begin{array}{*{20}l} {\sigma \left( {\frac{A}{{MR_{0}^{2} }}} \right) = \sqrt {\left( {\sigma C_{20} } \right)^{2} + \left( {\frac{{2\left( {2 - \gamma } \right)\sigma C_{22} }}{\gamma }} \right)^{2} + \left( {\frac{{24eC_{22} \sigma \theta }}{{\gamma^{2} \left( {3\theta - 6e} \right)^{2} }}} \right)^{2} ,} } \\ \end{array}$$

(8)

$$\begin{array}{*{20}l} {\sigma \left( {\frac{B}{{MR_{0}^{2} }}} \right) = \sqrt {\left( {\sigma C_{20} } \right)^{2} + \left( {\frac{{2\left( {2 + \gamma } \right)\sigma C_{22} }}{\gamma }} \right)^{2} + \left( {\frac{{24eC_{22} \sigma \theta }}{{\gamma^{2} \left( {3\theta - 6e} \right)^{2} }}} \right)^{2} ,} } \\ \end{array}$$

(9)

$$\begin{array}{*{20}l} {\sigma \left( {\frac{C}{{MR_{0}^{2} }}} \right) = \sqrt {\left( {\frac{{4\sigma C_{22} }}{\gamma }} \right)^{2} + \left( {\frac{{24eC_{22} \sigma \theta }}{{\gamma^{2} \left( {3\theta - 6e} \right)^{2} }}} \right)^{2} } .} \\ \end{array}$$

(10)

Figure 4a–c shows the errors of MOI for $A$, $B$, and $C$ normalized by $MR_{0}^{2}$, respectively. When creating these figures, we assume that $\sigma C_{20}$ and $\sigma C_{22}$ are the same and we used values of $C_{20}$, $C_{22}$, and $\gamma$ that are based on the homogenous model of Willner et al. (2014). From Fig. 4a–c, we can see that it is required to achieve a similar level of accuracy for both the libration amplitude and the gravity coefficients in order to realize a certain accuracy of MOI, otherwise the MOI error is dominated by either component that has the larger error. The required accuracy of the libration amplitude and the degree-2 gravity coefficient is 2% for A and B, and 3% for C.

Fig. 4

Error of moments of inertia (MOI) as a function of errors of degree-2 gravity coefficients and libration amplitude. Results for principal moments of inertia A, B, and C are shown from top (a) to bottom (c). Color bar indicates the MOI error in percent

After we defined the above-mentioned mission requirement using the simplified models in 2015, Le Maistre et al. (2019) presented more detailed and more realistic simulations using 6 different models of possible heterogeneous mass distributions inside Phobos. They computed geodetic observables based on different heterogeneity models by using discretized Phobos models with 500-m-length cubes that fit shape, volume, and total mass, and compared them with the observables based on homogenous mass distribution. Although our two-layer models are simple, they are complementary with the models proposed by Le Maistre et al. (2019) who did not explicitly investigate these two-layer models.

Le Maistre et al. (2019) show that the levels of departure of the MOI due to an interior heterogeneity are smaller than 15% of their respective homogeneous values, which is consistent with our results. Their evenly distributed buried ice (IED) models and disrupted and reaccreted (DAR-F1) models are comparable to our models shown in Fig. 2 (roch + ice) and Fig. 3 (rock + pore), respectively. The most striking departures from the homogenous case are found for models with mass excess (less porosity) or deficit (more porosity) beneath Stickney crater. They indicate that libration amplitude and low degree gravity coefficients at a precision better than 5% can allow to identify such heterogeneities, which can be detected if the MMX geodetic mission requirement is met.

Orbits and instruments for MMX geodetic study

Quasi-satellite orbit

The orbit of the spacecraft is fundamental information when observing a planetary body. The spatial resolution of the observation depends on the distance between the spacecraft and the target body. Observation plans are also made based on predicted orbit information. From the viewpoint of global observation, a polar orbit around Phobos is ideal. However, it is impossible for a spacecraft to naturally orbit around Phobos with a Keplerian motion, because Hill’s sphere of influence (SOI), that is the region around the body where the dynamics is dominated by its own gravity field, is small for Phobos (Zamaro and Biggs 2015). Therefore, a more stable quasi-satellite orbit (QSO) is proposed as the nominal orbit for MMX (Ikeda et al. 2017; Baresi et al. 2021). A QSO is a periodic orbit in the circular restricted three-body problem (CR3BP) and is a kind of distant retrograde orbit (DRO) (Lo and Paker 2004). In such a QSO, a spacecraft apparently goes around Phobos in a Phobos-fixed rotating frame, while it actually orbits around Mars with slightly different orbital elements from those of Phobos. A planar orbit is referred to as 2D-QSO (or just QSO for brevity), while one with non-zero out-of-plane component is referred to as 3D-QSO.

Figure 5 shows proposed planar QSOs with different sizes denoted by pericenter and apocenter distances. The QSOs are roughly divided into three categories, i.e., high (QSO-H, 100 × 200 km), mid (QSO-M, 50 × 100 km), and low (QSO-Ls < 50 km). The QSO-Ls are further divided into three sub-categories depending on their size (QSO-LA, LB, and LC). The planned QSO-H and -M orbits are nearly ellipses centered on Phobos with an axis ratio of about 1:2. The central distance of the spacecraft is minimum at sub-Mars and anti-sub-Mars points, and maximum at leading and trailing points. For operational safety reasons, the spacecraft will be inserted in a QSO-H first, then gradually lower its altitude. For the QSO-H configuration, the apparent orbital period in the Martian moon-centered Sun–Phobos synodic frame is about 35 days and we will be able to observe the dayside of Phobos for about 2 weeks. In order to provide a shape model during the early mission phase, the current plan is to conduct the main observation for shape modeling at the first QSO-H opportunity (see also Nakamura et al. 2021). Because the planned duration of the first QSO-H is 1 month, we investigate whether necessary observations can be completed within the 2-week dayside period described in "Shape modeling strategy at QSO-H" section.

Fig. 5

Five planar QSOs with different sizes plotted in Phobos-fixed rotating frame, projected on to the equatorial plane, and seen from the north pole. 1 day of orbits are shown. The spacecraft moves clockwise. Phobos’ shape is shown in red

Because the motion of a spacecraft in a 2D-QSO is almost confined to Phobos’ equatorial plane, only limited observations for high-latitude and polar regions can be made from 2D-QSO. However, this situation is greatly mitigated by adopting a 3D-QSO, thanks to its out-of-plane component. Currently, a 3D-QSO with the equivalent altitude of QSO-M is planned. Figure 6 shows an example of 3D-QSO with about 50-km out-of-plane excursion. The potential observation of polar regions from this 3D-QSO is discussed in "High-latitude and polar observations from 3D-QSO" section.

Fig. 6

Example of a 3D-QSO in Phobos-fixed rotating frame. A case with the in-plane size of 50 × 100 km and out-of-plane component up to Z = $\pm$50 km is depicted. The spacecraft position is shown in blue and Phobos’ shape is shown in red

Mission instruments for geodetic study

Three mission instruments and one bus instrument will contribute to generating geodetic products. These are the narrow-angle camera TENGOO, the wide-angle camera OROCHI (Kameda et al. 2021), the laser altimeter LIDAR (Senshu et al. 2021), and the narrow-angle camera CAM-T that is mainly used for spacecraft navigation purposes and has similar specifications as TENGOO. TENGOO has a field of view (FOV) of 0.82 × 1.09 deg and an instantaneous field of view (iFOV) of 5.8 $\mu$rad/pixel corresponding to 1.2 m spatial resolution at an altitude of 200 km. For comparison, HRSC onboard MEX has an iFOV of 9$\mu$rad/pixel (Oberst et al. 2008) corresponding to 1.8 m spatial resolution at an altitude of 200 km. Willner et al. (2014) used images from HRSC to reconstruct a shape model, of which spatial resolution varies between 4 and 80 m depending on flyby distances and channels. The high-resolution images obtained by TENGOO are to be used to update shape models, also for surface studies, landing site evaluation, etc. On the other hand, OROCHI has a wider field of view of 66 × 53 deg and an iFOV of about 400 $\mu$rad/pixel corresponding to 80 m spatial resolution at an altitude 200 km. If an OROCHI image is simultaneously taken with a TENGOO image, it provides important context information to the TENGOO image, i.e., which part of Phobos surface is imaged by TENGOO. Owing to the large field of view of OROCHI, Phobos can be framed within a single OROCHI image if the spacecraft altitude is larger than about 20 km. Therefore, OROCHI images will contribute to spacecraft navigation and orbit determination, i.e., spacecraft attitude information can be obtained from the limb shape and/or the position of landmarks. Because the field of view of TENGOO is narrow, a combination with CAM-T widens apparent field of view if the two cameras are configured with appropriate alignment offset. Such a combination will contribute to reducing the time duration required to cover a certain area of Phobos. This possibility is discussed in "Spatial and temporal coverage analysis" section.

LIDAR measures the line-of-sight distance between the MMX spacecraft and Phobos’ surface by measuring the round-trip travel time of the laser pulse. Its designed measurement range is from 100 m to 100 km assuming a geometric albedo of 0.05. Note that Simonelli et al. (1998) reported a globally averaged albedo value of 0.071$\pm$0.012. The divergence angle of the transmitted beam is 0.5 mrad, which corresponds to a footprint size of 50 m at an altitude of 100 km. The distance information provided by LIDAR will contribute to improving ephemeris, gravity field, rotational state, and shape of Phobos, as well as the spacecraft orbit. LIDAR range data provide a scale to image-based shape models and cumulative along-track LIDAR footprint (or laser bounce point) positions themselves will provide us high-fidelity geometric height information of Phobos topography. The accuracy of LIDAR measurement almost linearly varies with altitude from $\pm$2 m at an altitude of 100 m to $\pm$22 m at an altitude of 100 km (Senshu et al. 2021). Because the dominant part of the altitude-dependent error is a bias component, the relative along-track topographic profiles can be obtained with an accuracy better than 22 m even at an altitude of 100 km. The LIDAR measurements will contribute to verifying shape models in terms of topographic features whose horizontal scale is larger than the size of the LIDAR footprint. In addition to data from the onboard instruments, we will also use ground-based radio-tracking data (range and range rate) for geodetic study.

Although it is not directly related to geodetic observation, the mass spectrum analyzer (MSA) (Yokota et al. 2021) may contribute to the interpretation of the geodetic observations in terms of internal structure. Fanale and Salvail (1989, 1990) theoretically predicted that the current H₂O loss rate can be up to 3 g/s ($\sim$10²³ molecules/s) depending on porosity and pore size. Such sublimation will generate a H₂O gas torus along Phobos’ orbit (Mura et al. 2002). MSA can detect such water molecules coming from Phobos’ interior if the loss rate is larger than 10²² molecules/s. Fanale and Salvail (1889, 1990) also calculated the latitude-dependent water ice depth that also depends on porosity and pore size. If MSA detected water molecules from Phobos’ interior, this information would help us construct interior models of the water ice distribution based on the MSA observation, geodetic observations, and the theoretical studies.

Strategies for observations and data analysis

Shape modeling strategy at QSO-H

Because landing on Phobos is scheduled around the middle of about 3-year staying period around Mars (Nakamura et al. 2021), shape models are to be generated from images collected by the MMX spacecraft from a QSO during early mission phases. This section describes the currently planned strategy for MMX shape modeling by images from QSO-H. In contrast to the current situation of having a limited number of Phobos images with inhomogeneous spatial resolution (Karachevtseva et al. 2014; Witasse et al. 2014), MMX will provide a larger number of images with more homogenous quality and resolution. An existing shape model will be used as initial base model to which our first update will be made with telescopic camera images to be acquired during the terminal rendezvous phase from a few thousand kilometers distance, with entire or part of Phobos being framed within the field of view of the camera. Due to Phobos’ synchronous rotation, only certain hemisphere can be imaged, but these new data would nevertheless contribute to update the initial model. The following two subsections describe preliminary coverage analyses for topography reconstruction from images taken from QSO-H.

Spatial coverage analysis

Because spacecraft nadir nominally directs to the equator during 2D-QSO periods, we need to slew the spacecraft attitude in latitudinal direction to observe non-equatorial regions. We will first focus on covering equator to mid-latitude region. Note that possible landing site candidates are to be selected within $\pm$30° latitudinal band. Figure 7a shows an example of multiple TENGOO field of view to cover a $\pm$45° latitudinal band which is calculated based on a simulated QSO-H trajectory, the TENGOO field of view parameters, and planetographic latitude calculated with a flattening of 0.3. Note that we will make the long side of the rectangle field of view aligned with cross-track (or latitudinal) direction in 2D-QSO operations. The $\pm$45° latitudinal band is covered by 5 to 9 TENGOO images for a certain longitudinal section depending on the longitude (or spacecraft altitude). In total, 508 images from 66 longitudinal sections are necessary to cover the latitudinal band with a 15% of FOV overlap between adjacent images. When similar observations are made from QSO-M, these numbers increase to 8–20 TENGOO images per longitudinal section and 2459 images from 154 sections (not shown).

Fig. 7

Simulated field of view of the telescopic camera TENGOO projected on Phobos’ surface. Top (a) coverage achieved by TENGOO and bottom (b) coverage achieved by an imaginary camera that mimics a combination of TENGOO and CAM-T with effective longitudinal field of view almost doubled. Observations from equatorial QSO-H with an off-nadir pointing to cover a $\pm 45^\circ$ latitudinal band are planned

As will be shown in "Spatial and temporal coverage analysis" section, the first dayside observation opportunity from QSO-H will be too short to cover enough images for shape modeling. One idea to overcome this situation is to spatially complement the insufficiency by making use of CAM-T together with TENGOO. If the fields of view of the two cameras are configured as indicated in Fig. 8, it will result in almost doubled viewing width for both the along-track (longitudinal) and cross-track (latitudinal) directions when scanned in cross-track and along-track directions, respectively. For the current case of taking images of non-equatorial regions from 2D-QSO by slewing the spacecraft attitude in cross-track direction, the camera alignment offset in the along-track direction takes effect in extending longitudinal field of view. Figure 7b shows the result of a simplified coverage analysis where the viewing width is set to1.52° instead of 0.82° assuming a FOV overlap of 15%. The number of longitudinal sections needed is reduced from 66 to 33 and the total number of images needed to cover the regions of interest is reduced from 508 to 262. In fact, the total number of images does not change so much because TENGOO and CAM-T individually provide images, but the wider area covered by a single favorable observation opportunity will make it possible to complete the imaging operations required for shape modeling within the first dayside period of QSO-H (see next subsection).

Fig. 8

Proposed mutual alignment between TENGOO and CAM-T

Spatial and temporal coverage analysis

This subsection focuses on whether observations from QSO-H that are necessary for shape modeling are completed within the first dayside period of about 2 weeks. The shape model construction technique of SPC is one of the candidates to be applied to MMX. This technique has been used for numerous celestial bodies, e.g., 25143 Itokawa (Gaskell et al. 2008), 162173 Ryugu (Watanabe et al. 2019), 1 Ceres (Park et al. 2019), Phobos and Deimos, (Ernst et al., submitted to Earth Planets and Space) etc. Here we assume that this technique is applied for MMX shape modeling. Because SPC requires variations in both emission (viewing) and incidence (illumination) angles, we plan to observe a certain area on Phobos surface at least seven times with the following conditions: (1) three times from east, west, and the meridian of the target area at local time around 9.5 h; (2) three times from east, west and the meridian directions at local time around 14.5 h, and (3) once from the meridian around noon (12.5 h).

We search for timings of favorable observations that meet the above conditions within certain tolerance, i.e., we focus on a surface point on the equator at the target’s local time and calculate the emission angle at a certain epoch. On the other hand, the whole Phobos is split into 60 sectors each covering a longitude interval of 6° that corresponds to the longitude width covered by TENGOO at the lowest altitude of QSO-H. The temporal resolution of the analysis is set to 10 min, considering the time required for the cross-track attitude slewing plus margin. The spacecraft position at a certain epoch is judged as “favorable” when the emission angle is 20$\pm$5 degrees at the equator. With this condition met, the longitudinal sector, on which the equatorial point at target local time falls, is marked as “observed”. The spacecraft cannot observe Phobos all the time because the spacecraft-fixed high-gain antenna is supposed to point to Earth during periods reserved for bus operations including telecommanding, spacecraft tracking and science/engineering telemetry downlink. The required time is estimated to be 8 h/day, which is excluded from our analysis as “unobservable” period that is set from 0 to 8 UT considering the visibility condition of the spacecraft from a tracking station in Japan. Plotted in Fig. 9 are the timing for the seven “favorable” observation opportunities and “observed” longitude, which shows that the observation opportunities are temporally well separated over the dayside period. However, there are many longitudinal sectors that are not observed seven times during the first QSO-H dayside period as shown in Fig. 10a. The situation is more severe for longitudes around 0 and 180° (Mars and anti-Mars sides) where the spacecraft’s altitude becomes minimum and the resultant surface area covered by TENGOO’s FOV becomes the smallest. The spacecraft velocity relative to Phobos’ surface is also maximum at these longitudes. Next, we analyzed the “TENGOO plus CAM-T” case where the width of the longitudinal sector is set to 12° and the result is shown in Fig. 10b. In this case, almost all of the longitudinal sectors are observed with the seven measurement opportunities, which means that the necessary observations can be mostly completed within the first dayside period by combining TENGOO with CAM-T. We plan to fill the remaining observation gaps in the subsequent QSO-M period (see Fig. 1 of Nakamura et al. 2021).

Fig. 9

Favorable observation timings for 7 different conditions. Pluses: looking at local time (LT) 14.5 h from east, crosses: looking at LT = 14.5 h from the meridian, asterisks: looking at LT = 14.5 h from west, open boxes: looking at LT = 12.5 h from the meridian, filled boxes: looking at LT = 9.5 h from east, open circles: LT = 9.5 h from the meridian, filled circles: looking at LT = 9.5 h from west. The viewing angle (emission angle) from east and west is set to 20° at the equator. Small dots indicate the nadir longitude of MMX spacecraft every 10 min. Data during the bus operation period are excluded from the plot

Fig. 10

Number of observation conditions (up to 7) realized for each longitudinal sector: a TENGOO only case and b TENGOO plus CAM-T case. The latter case assumes that the longitudinal FOV is doubled

The above analysis is based on a simplified spheroidal shape of Phobos. The real topography of Phobos shows large deviations from the spheroid, e.g., inside Stickney, a large crater of 9 km across. A more detailed analysis is necessary for such regions because incidence and emission angles show complex variations due to the topographic slopes.

Images with higher spatial resolution are to be obtained from QSOs with lower altitude than QSO-H. The survey for landing site candidates is planned from QSO-M, and more detailed observations will be made from QSO-Ls for further down-selected candidates. Local digital terrain models will be developed by SPC and/or stereo photogrammetry using TENGOO/CAM-T images from these low-altitude QSOs.

High-latitude and polar observations from 3D-QSO

A 3D-QSO provides us with precious opportunities to observe high-latitude and polar regions of Phobos. The observability depends on the amount of the out-of-plane excursion that correlates with delta-V amount required to transfer from 2D-QSO. Currently an excursion from planar QSO-M up to about 50 km is proposed (Fig. 6). Figure 11a shows the nadir-pointing ground tracks with the latitude being planetographic (flattening = 0.3). For this particular 3D-QSO, the ground track repeat cycle is about 9 days. On the other hand, the nadir illumination conditions vary with a period of about 1 week and it is found that, if the observation is restricted to the dayside region only, the corresponding apparent repeat cycle is about 1 month.

Fig. 11

Top (a) nadir-pointing spacecraft ground tracks for the 3D-QSO shown in Fig. 6. Bottom (b) ground tracks making polar observations possible by formulating the slewing angle as a function of the Z-component of spacecraft position. The orbits are plotted for a span of 10 days

Figure 11a indicates that an off-nadir pointing is necessary to look at the poles. The required slewing angle depends on the spacecraft position. Figure 11b shows an example of a ground track that makes polar observations possible by formulating the slewing angle as a function of Z-component of spacecraft position (or the out-of-plane component). The illumination condition of the polar regions varies with a beta angle that has seasonal variations with an amplitude of about 25°, where the beta angle is the angle between the orbital plane of Phobos and the vector to the Sun or the latitude of the sub-solar point on Phobos. The MMX spacecraft will be inserted into QSO-LB or LC when the south polar region is illuminated with a beta angle about − 25° around May 2026 and landing operations is planned around April 2027 when the north polar region is illuminated with a beta angle of 25°. The envisioned time frame of the 3D-QSO is around October 2026 when the beta angle is close to zero (Nakamura et al. 2021). The beta angle would change by about 6° if the timeline shifted by a month. Such a low absolute value of the beta angle is suitable to get a first overview of both poles. However, the incidence angle becomes large in the polar regions, which makes topography reconstruction from images difficult. Although its spatial resolution is lower than that of images, LIDAR data will help reconstruct the topography of polar regions. The maximum slant range is about 120 km that exceeds the designed LIDAR range capability of 100 km, but the designed capability is based on a geometric albedo of 0.05 that is smaller than the global average of 0.07. If the polar albedo is equal to or larger than the global average, the distance attenuation of reflected laser power can be compensated for, making LIDAR observation possible, albeit not perfect due to yet large incidence angle of the LIDAR measurements.

Strategies for ephemeris, gravity field, and rotation parameter of Phobos

Because the orbit plane of the MMX spacecraft is close to that of Phobos having a small inclination with respect to Mars equatorial plane, the spacecraft velocity change would not be perpendicular to the line-of-sight direction to the Earth. It is expected that spacecraft orbit around Mars will be well determined from Earth-based range and Doppler observations with good sensitivity owing to the above geometry. The Phobos ephemeris can be improved by combining the spacecraft position in the inertial frame determined by the range and Doppler data with target-relative LIDAR and image data. Most significant error in the mean motion (Jacobson and Laney 2014) is expected to be reduced. Phobos ephemeris improvement will be enabled by operating LIDAR from QSO-M and -L with even better performances if radio-tracking of the spacecraft (range and Doppler measurements) using medium-gain antenna is carried out simultaneously with LIDAR ranging. Yang et al. (2019) estimated errors of degree-2 gravity coefficients introduced by error in Phobos ephemeris, which was no larger than 0.52% based on an assumed shift of 200 m. If true ephemeris error was less than about 1 km, the corresponding gravity coefficients error would be smaller than the required accuracy.

Because the gravity field perturbs the spacecraft orbit, a precise orbit determination improves the estimation of the gravity coefficients. From the viewpoint of gravity field recovery, a smaller size of QSO or lower spacecraft altitude is preferable because a larger sensitivity to the gravity field is gained by shorter distance between the spacecraft and the target. In this regard, QSO-H provides only weak sensitivities to gravity coefficients because the size of QSO-H is larger than 9 times of Phobos mean radius ($\sim$11 km). Spacecraft tracking at QSO-M will allow us to better retrieve $C_{22}$, and 3D-QSO at this altitude will contribute to solving for zonal $C_{20}$ with smaller correlation with $C_{22}$. Further observation at QSO-Ls, in particular the planned lowest one QSO-LC of 20 × 27 km in size with a minimum altitude of about 7 km, provides a better opportunity for estimating $C_{22}$ and possibly other low-degree sectorial coefficients. The best situation for the gravity field recovery would be given by a 3D-QSO at such a low altitude, but the feasibility depends on the stability of the orbit and the safety of the spacecraft.

In addition to the conventional Earth-based 2-way Doppler and range observation, LIDAR data and optical landmark tracking data are to be incorporated into the orbit determination process to improve its accuracy. The LIDAR observation provides target-relative line-of-sight distance to Phobos, and the optical observation provides pixel and line location of a landmark in the image which has sensitivity in the direction perpendicular to the line-of-sight direction. The two observations play complementary roles in the orbit determination process (Yamamoto et al. 2020).

The acceleration on the spacecraft caused by the gravity coefficients is calculated in the body-fixed coordinate system. If the axes of the body-fixed frame are defined by the principal axes of the body, $C_{21}$, $S_{21}$ and $S_{22}$ disappear. A precise shape model is needed to translate the target-relative LIDAR measurements to the Phobos center of mass (COM) that is necessary for the spacecraft orbit determination and Phobos ephemeris determination. The offset between center of mass and center of figure (COF) will be first calculated from the image-based shape model by assuming constant density, based on which the orbit and gravity analysis will be conducted. The resultant degree one gravity coefficient will be used to redefine the COM–COF offset for subsequent iterative analyses.

The main error sources for a SPC shape model are the uncertainties in the camera position (or spacecraft position) and camera pointing direction, which are difficult to separate because the image data alone are not sensitive to radial position of the camera (Park et al. 2019). Errors in rotation parameters (e.g., libration amplitude) will also degrade the shape modeling accuracy. If the orbit determination process can estimate the spacecraft state, gravity field, rotation parameters, and camera pointing directions, the shape modeling results will be improved by these information, which will in turn update the orbit determination results by iterative procedures. Such geodetically derived information will contribute to selecting potential landing sites as well as to the planning of the landing operation itself.

Summary

This paper described the accuracy of libration amplitude and degree-2 gravity coefficients that is necessary to meet a part of MMX mission requirements. We studied how the MOI changes due to an inhomogeneous distribution of water ice or porosity with respect to the MOI for a homogenous interior of Phobos. The calculation is based on a simple two-layer model having the similar boundary figure to the surface, which is consistent with the currently observed libration amplitude. In such a configuration, the maximum deviation of the MOI is smaller than 17%. In order to detect the presence of inhomogeneity in the density structure due to localization of water ice of 10% of Phobos’ mass, it is required to achieve at least 3% of MOI accuracy. To this end, the required accuracy for the libration amplitude and the degree-2 gravity coefficients are 2–3%.

The data from three cameras of TENGOO, OROCHI, CAM-T and laser altimeter LIDAR, together with Earth-based radio-tracking data, will contribute to improving our knowledge of ephemeris, gravity field, shape, and rotation parameter of Phobos. We presented the preliminary imaging strategy for creating a base shape model of the equatorial to mid-latitude region of Phobos that is required for landing site selection process at the early stage of the mission. The accuracies of the geodetic products are expected to be improved by an iterative process of orbit determination and shape modeling. Such products will contribute to achieving the MMX mission goals by providing basic and important information for interpreting data from MMX mission instruments and for selecting the landing site.