Simulation of gravity field estimation of Phobos for Martian Moon eXploration (MMX)

We study the dynamic orbit about Phobos of the Martian Moon eXploration (MMX) spacecraft and simulate gravity field estimation using Doppler, image landmarks, and LIght Detection And Ranging (LIDAR) data based on a mission plan of MMX and investigate whether the differences in the internal density structure of Phobos can be detected through this mission. Degree 2 Stokes coefficients of the Phobos gravity field, C 20 and C 22 , which are necessary to obtain the moment of inertia by combining with the libration amplitude, could be determined with an accuracy within the order of 0.1% using datasets from two-dimensional (2D) quasi-satellite orbits (QSOs). If observations from 3D-QSO could be realized at low altitude, coefficients up to degree and order 5 could be estimated, which could be used to detect regional density anomalies. Moreover, the observation data from the ascending trajectory after landing on the Phobos surface can be used to detect local density anomalies around the landing site, which could help interpret the origin of the samples and understand the internal structure of Phobos.


Graphical Abstract 1 Introduction
Mars is located at the edge of the inner solar system, and liquid water is believed to have been widespread on its surface in the past.Thus, studying the formation and evolution of Mars has gained attention.However, traces on the Martian surface for formation and evolutionary processes are not evident and are insufficient to acquire information about the formation and evolution of Mars.Contrastingly, the surfaces of the Martian moons are comparably pristine and are likely to have preserved the past records.Thus, studying the Martian moons may provide important information about the evolutionary histories not only of themselves, but also of Mars.Following the successful Hayabusa and Hayabusa2 asteroid missions, the Martian Moon eXploration (MMX) mission is the third sample return mission conducted by the Japan Aerospace Exploration Agency (Kuramoto et al. 2022) scheduled to launch in 2026 and return to Earth in 2031.The target bodies of the MMX mission are the two Martian moons, Phobos and Deimos.
The origin of Phobos and Deimos is still controversial.There are two major hypotheses for their origin: capture of asteroids and in situ formation by a giant impact on Mars.The capture hypothesis proposes that primordial celestial bodies, which were formed outside the snow line, moved inward and subsequently, became satellites as Mars captured them.The low surface reflectance of Phobos and Deimos and their reddish spectra without clear silicate absorption are similar to D-type asteroids, which are considered a strong argument to support the capture hypothesis (Burns 1978;Murchie et al. 1991;Pajola et al. 2013).If the hypothesis is correct, the MMX mission will provide significant insights into how water and other volatiles were transported from the outer Solar System to terrestrial planets.Conversely, the currently observed orbital eccentricities and inclinations of the Martian moons with respect to the Martian equatorial plane are small.Moreover, the orbit of Deimos is outside the synchronous orbit of Mars, while that of Phobos is inside of the Mars orbit.This can be explained more naturally by considering the in situ formation of the satellites by a giant impact on Mars (Ida et al. 1997;Citron et al. 2015;Hyodo et al. 2017;Canup and Salmon 2018).An extended giant impact hypothesis has also been proposed, which suggests that Phobos underwent multiple tidal destructions since its formation (Hesselbrock and Minton 2017).
The Martian moons have rarely been directly explored so far, and most previous observations were from Mars explorers, which occasionally conducted fast flybys or which moved in low-Mars orbits and conducted observations from large distances (which is particularly true for Deimos).However, MMX will conduct comprehensive remote sensing observations in the vicinity of Phobos, along with sampling and rover observations on its surface, while flyby observations have been planned for Deimos.Further, continuous observations of the Martian atmosphere have been planned.Using the obtained precise datasets, elucidating the evolution of the Mars system and Solar system is one of the main scientific objectives of the mission (Kuramoto et al. 2022).
In this study, we focus on estimating the gravity field of Phobos through the MMX mission.The celestial parameters of Phobos are shown in Table 1.The gravity field of a celestial body provides important information to infer its internal mass distribution, which reflects the origin and evolutionary process of the body.Using radio-tracking and image datasets acquired during Phobos flyby of Mars Express, the mass value (Rosenblatt et al. 2008) and shape model (Willner et al. 2010(Willner et al. , 2014) ) were significantly improved from previously acquired values and the resulting bulk density of Phobos was estimated as 1850 ± 70 kg/ m 3 (Rosenblatt 2011).However, this value is smaller than that of most carbonaceous meteorites, and thus, Phobos is expected to contain light elements, high porosity (void), and/or water ice in its interior (Willner 2009).Based on the estimated density and hypotheses about the origin, Le Maistre et al. (2019) considered several possible large-scale internal density structure models of Phobos, and discussed the precision of the gravity field, rotation, and tidal parameters required to distinguish these differences.Furthermore, Matsumoto et al. (2021) aimed to distinguish the differences in the radial-direction internal density structure models of Phobos using the moment of inertia (MOI) values and showed the required precision of the libration amplitude and degree-2 gravity field coefficients from which MOI is calculated.
A concerning issue is whether the parameters required to distinguish such large-scale internal structures can be estimated with sufficient precision in the actual MMX mission.Hence, in this study, realistic spacecraft orbit and observation datasets, based on the MMX mission plan, are prepared and used to simulate the estimation of the spacecraft orbit and Phobos gravity field in the MMX mission.
The remainder of the paper is organized as follows: Sect. 2 describes the development of alternative models for simulating the internal structure of Phobos.Realistic but simulated spacecraft orbit and observation data are also generated based on the observation plan of the MMX mission.These datasets are used to simulate orbit and gravity field estimations described in Sect.3. The simulation results are presented in Sect. 4. In Sect.5, whether the differences in the various models of the internal density structure of Phobos are detectable are discussed, and procedures to achieve the interior structure models of Phobos more efficiently are recommended in the frame of the mission plan.Finally, Sect.6 provides the conclusions.

Internal density structure model of Phobos
Based on the hypotheses of the formation, various simple internal density structure models of Phobos were proposed for simulating internal structure estimation based on geodetic parameters (Le Maistre et al. 2019;Dmitrovskii et al. 2022).Among these, typical models, for example, are homogeneous density model, which assumes a rubble pile body (Murchie et al. 1991), a model with radially layered structures with different density substances (e.g., rubble, monoliths, ice, or voids for each layer) (Britt et al. 2002;Fanale and Salvail 1990), and a model with density anomaly below the largest crater Stickney due to soil compaction during crater formation (Britt et al. 2002;Richardson et al. 2002).Referring to these models, we develop some Phobos models, as described below, with different internal density structures for our simulation.
We assume a fully uncompensated Phobos model, and do not consider mass compensation.The mechanism by which celestial bodies support their topography depends on their size as well as their internal structure.On a geological timescale, rocky and icy bodies with a radius exceeding about 200 km typically exhibit shapes that are close to the shape of hydrostatic equilibrium but can maintain some non-hydrostatic topography.The mechanism that associates gravity with topography involves the elastic properties of the lithosphere, the viscous stresses from mantle convection, and particularly over long-time scales, Airy and Pratt isostasy with respect to the loading of large-spatial-scale topography or the loading created at the time of formation of the lithosphere (Hemingway and Matsuyama 2017).However, for minor bodies with radii smaller than 200 km, the shapes often significantly deviate from the ones in hydrostatic equilibrium, and it is believed that the mechanism supporting the surface topography is predominantly governed by frictional strength.Thus, for minor bodies, uncompensated topography is usually a good assumption (Ermakov et al. 2018).

Homogeneous density model
In this study, the shape model created by Willner et al. (2014), consisting of 274,874 facets and 137,439 vertices, Table 1 Mean celestial parameters of Phobos at J2000 (TDB) referred to the local Laplace plane derived from MAR097 ephemeris (Jacobson and Lainey 2014) a: semi-major axis; e: eccentricity; ω : argument of periapsis; M a : mean anomaly; i 0 : inclination with respect to the Laplace plane; : longitude of the ascending node on the Laplace plane; P: sidereal period; P ω : argument of periapsis precession period; P : longitude of the ascending node precession period; α : right ascension of the Laplace-plane pole position with respect to ICRF; δ : declination of the Laplace-plane pole position with respect to ICRF; ǫ : angle between the planet's equatorial and Laplace planes is used for the simulation.To develop the homogeneous density model, we first approximate the shape model as a sum of equal-sized small cubes (Willner 2009), with the size of the cubes set to 50 × 50 × 50 m.The total volume of this approximate shape model is 5748.00km 3 .The density of each small cube is set to 1857.82 kg/m 3 , which is obtained by dividing the total mass 1.068 × 10 7 kg, which is based on the standard gravitational parameter (GM) value determined by Andert et al. (2010), by the total volume.

Layered density model
The density distribution of the inhomogeneous Phobos is developed by changing the density of each small cube under the constraint of volume and mass conservation.For layered density models, in principle, numerous models can be developed by changing the number of layers, density of each layer, and thickness and depth of each density boundary.However, this study does not aim to create the most plausible Phobos internal structure model, but to investigate the precision of the estimated parameters required to distinguish the different internal density structures.Thus, we create three simple radially non-uniform two-layer density structure models, which are described in Table 2.We assume that the density boundary is 80% of the radial distance from the center to each surface facet that composes the shape model, and that the boundary shape is similar to the surface profile.
The water contents of the three models are assumed to be 10, 20 and 30% of the total volume, respectively, and the outer layer is assumed to comprise water and rock, while the inner layer is assumed to be composed of only rock.Moreover, we assume that the outer layer is a homogeneous layer consisting of a mixture of ice and rock to conduct simple simulation assuming a two-layered density structure.In reality, however, it is more plausible that there is another thin surface layer of rock on top of the ice (or mixed ice and rock) layer (Fanale and Salvail 1990), as ice sublimates quickly on the surface.

Regional density anomaly model
A model with density anomaly beneath the Stickney crater is developed under the constraint of total mass and volume and the amplitude of physical libration in longitude.The libration amplitude of Phobos was obtained with a formal error of 3% from previous studies (Burmeister et al. 2018;Lainey et al. 2021), but the values differed between the studies.Thus, constraining the internal structure of Phobos by the current knowledge of the libration amplitude is insufficient and the value is to be refined in the MMX mission (Matsumoto et al. 2021).Burmeister et al. (2018) obtained the libration amplitude of Phobos (− 1.14° ± 0.03°) by bundle block adjustment using the images from the Mars Express and Viking missions.This value implied that the libration amplitude was mostly caused by the highly irregular shape of Phobos, and that there were no significant regional density anomalies in its interior.Conversely, the libration amplitude (− 1.09° ± 0.01°) (Lainey et al. 2021), which was derived by analyzing the Phobos and Deimos ephemerides data, allowed for the presence of regional density anomalies.Herein, we set the libration value to − 1.08°, assuming a 2-sigma error from the value derived by Burmeister et al. (2018), and use it to create the local density anomaly model of Phobos.
If Phobos consists of the sum of small cubic masses with two different densities, ρ 1 and ρ 2 , the normalized moment of inertia (MOI) of Phobos about each of the three principal axes can be expressed as follows: where a 1 , a 2 , b 1 , b 2 , c 1 , and c 2 are given as: where A , B , and C are the normalized MOIs of Pho- bos about the x, y, and z axes, respectively.x i , y i , z i and V are the center coordinates and volume of each small (1) cubic mass element, respectively, and N 1 and N 2 are number of small cubic masses with ρ 1 and ρ 2 densities, respectively.Moreover, M is the total mass of Phobos that is equal to ρ av (N 1 + N 2 )�V , where ρ av is the average density of Phobos, and R 0 is the mean radius of Phobos.Using Eq. ( 1), the dynamic flattening of Phobos, γ = B−A C , can be written as: The conservation equation of the total mass of Phobos is: Using Eqs. ( 2) and (3), the densities ρ 1 and ρ 2 can be derived as follows: Additionally, the libration amplitude θ is related to γ as: where e is the orbit eccentricity of Phobos, which is set to 0.01511 (Jacobson and Lainey 2014) in this study.
According to Eqs. ( 4) and ( 5), if the spatial distribution of the density inside Phobos is specified, ρ 1 and ρ 2 values can be determined under the constraints of the total mass, volume, and libration angle.In this study, as a model with density anomaly beneath the Stickney crater, we assume that there is a density anomaly region of the horizontal radius 20° centered at (longitude, latitude) = (50°W, 0°N).The depth of the density boundary is assumed 70% of the radius from the center to the surface.In this case, the density is 2827.93 and 1836.15kg/m 3 for the density anomaly region and other regions, respectively.Figure 1 shows the gravity anomaly map calculated from the developed regional density anomaly model.

Development of spherical harmonic coefficients of gravity field
The fully normalized spherical harmonic coefficients of the gravity field (Stokes coefficients) corresponding to each internal density structure model are obtained using the method of Werner (1997) from each of the Phobos models comprising small cubes.The coefficients up to degree and order 180 are used in this study.The origin of the reference frame is defined as the origin of the shape model, and the rotational axis is aligned with its z-axis.The shape model is created based on imaging observations from past missions, and the origin corresponds (3) to the observed center of mass.On the other hand, the center of figure of this shape model, which is calculated as center of mass with constant density assumption (centroid), does not coincide with the origin of this shape model and has an offset of approximately 125 m.With our current knowledge, it remains unclear whether this difference is due to errors in shape modeling or is a result of density inhomogeneities.Therefore, we do not make any adjustments to correct the difference between the center of mass and center of figure and left it as is.Thus, the degree 1 terms of the generated spherical harmonics for each model are not zero.

Simulated orbit of MMX spacecraft
The center of orbit integration is set to the Phobos center.
For the orbit integration, we consider the gravity field force of Phobos, solar radiation pressure force, and point mass forces of Sun, Mercury, Venus, Earth, Mars, and Jupiter.For Mars, in addition to the point mass force, the perturbing force due to the higher-degree gravity field is employed using the coefficients up to degree and order 20 of the Goddard Mars Model 3 (Genova et al. 2016).
For the gravity field model of Phobos, Stokes coefficients up to degree and order 180 developed using the procedure mentioned in Sect.2.1 are used.At lower altitudes of the ascending trajectory (Sect.4.2), when the spacecraft is positioned inside the Brillouin sphere of the spherical harmonic expansion, we do not use the coefficients to calculate the gravity force but calculate it as the sum of the gravity force from each small cubic mass from the Phobos density structure model.The size of each cubic mass is 50 × 50 × 50 m, while the maximum degree/ order of Stokes coefficients used is 180.Hence, in principle, the gravity acceleration derived by the summation method includes shorter wavelength components than those derived from the spherical harmonic methods.This might cause a discontinuity in the gravity acceleration value at the switching altitude of the two calculation methods.Therefore, we compare the gravity accelerations calculated by the two methods at the Brillouin radius using the homogeneous density Phobos model.The results shows that the difference between the two was 0.003%, confirming that the orbital discontinuity is almost negligible.
MAR097 (Jacobson and Lainey 2014) is used for the Mars and Phobos ephemerides, and DE430 (Folkner et al. 2014) is used for the Sun and other planets.For the rotation parameters of Phobos, Mars, and Earth, the models described in the 2009 IAU Working Group report are used (Archinal et al. 2011).A simple model including two solar panels and one box is used as the MMX spacecraft model to calculate solar radiation pressure force.The reflectivities of the solar panel and box are set to 0.010 and 0.375, respectively.It is assumed that the solar panels are always oriented towards the Sun during observations.
The MMX spacecraft is scheduled to be launched in 2026 and will arrive in the Mars orbit in 2027.
The spacecraft will initially co-orbit with Phobos, which would have moved to a phasing orbit, and finally enter a quasi-satellite orbit (QSO) from which Phobos would be observed.QSO is an orbit around Mars, and also around Phobos when viewed in the Phobos-fixed frame.The spacecraft is scheduled to stay around Mars for three years, with operations planned for the following five scientific observation phases.Phase 1: initial checkout of instruments and first observations of Phobos, Phase 2: Phobos proximity observation to select landing sites, Phase 3: two landings on Phobos and delivery of the rover to Phobos, Phase 4: scientific observations of Phobos and Mars, Phase 5: scientific observations of Deimos and Mars.(Nakamura et al. 2021).In Phase 1 and 2, the QSO operation will begin at a high altitude for the safety of spacecraft operations and will gradually move to lower altitudes while updating the physical parameters necessary for the operation.In addition to the observations from the QSOs whose orbital plane is almost in the equatorial (2-dimensional, 2D) plane of Phobos, observations from the 3D-QSOs with a large orbital inclination angle are also proposed as one of the possible orbits, to cover high-latitude areas.In this study, based on these orbit plans, some QSOs with different altitudes and latitudinal coverages, i.e., QSO-H, QSO-LA, 3D-QSO-M, and 3D-QSO-LA are generated and used for the simulations, although the 3D-QSOs are not nominal as of this writing.In a Phobos-fixed frame, the spatial dimensions for the QSO-H, QSO-LA, 3D-QSO-M, and 3D-QSO-LA are ± 200 × ± 100 × ± 4, ± 50 × ± 30 × ± 0.4, ± 50 × ± 100 × ± 40, and ± 50 × ± 30 × ± 10 km, respectively (Fig. 2).In the MMX mission, the spacecraft is also planned to temporarily descend from the QSO-LA altitude to collect surface samples or drop a rover on the Phobos surface, after which it will ascend and return to the QSO-M altitude.During the descent, the orbit will be frequently controlled owing to the strict requirements for landing or rover insertion sites, resulting in significant perturbations of the spacecraft trajectory.Therefore, the descending trajectory is not suitable for orbit analysis to detect the gravity anomaly of Phobos.Conversely, the ascending trajectory will be ballistic with no orbit control, which is a good opportunity to observe the Phobos gravity field at low altitudes.Therefore, in addition to QSOs, this study also generates an ascending spacecraft trajectory and investigate the possibility of detecting the Phobos density anomaly from this trajectory.We assume an initial position 100 m above the Phobos surface at (longitude, latitude) = (0°E, 0°N) and an initial velocity of 10 m/s upward with respect to the surface, and generate an ascending ballistic trajectory until the spacecraft reaches an altitude of 50 km.

Simulated observation data
Based on the observation plan of the MMX mission, we assume that the following observation datasets are available for estimating the orbit and gravity field: Doppler data derived from ground radio-tracking observations, range data derived from on-board laser altimeter observations, and landmark coordinate data on the image plane acquired through observations captured by on-board cameras.The simulated observation datasets are created as follows.
We assumed that the Doppler observation data are acquired by radio tracking of the MMX spacecraft at elevation angles of 20° or more.In actual data processing, the obtained Doppler counts are converted to timeaveraged range rates after acquisition (Montenbruck and Gill 2000); however, in this study, for the simplicity of the simulation, instantaneous range rate data are acquired instead of Doppler data.The range rate data interval is set to 60 s, assuming that the obtained Doppler counts are integrated every 60 s.
The observation range data between the MMX spacecraft to the Phobos surface are acquired by on-board laser altimeter, LIght Detection And Ranging (LIDAR) system (Senshu et al. 2021).The range data used in the simulation are created from the spacecraft position and attitude, and laser footprint coordinate on the Phobos surface at each time.To simplify the simulation, we assume that the attitude each time is set so that the laser line-of-sight vector is oriented toward the Phobos center.The laser footprint coordinate is calculated using the SPICE toolkit (Acton et al. 2018), and the data interval is set to 1 s.
The landmark coordinate data on the images used in this study are created as follows: initially, 3D landmark coordinates in the Phobos-fixed frame are defined on the Phobos shape model.Later, 410 globally distributed landmarks are defined and used.However, there may be cases that none of these landmarks are included in the image captured by a telescopic camera at low altitudes.In such cases, we use an increased number of the defined landmarks (40,924) for the simulation, by expecting that low-altitude image observations can provide a larger number of landmarks on an improved resolution shape model.Subsequently, we calculate the 2D pixel coordinates of these landmarks at each time taken by two onboard cameras, namely, a wide-angle camera (Optical RadiOmeter composed of CHromatic Imagers, ORO-CHI) and a telescopic camera (TElescopic Nadir imager for GeOmOrphology, TENGOO), and use them for the simulation.The camera parameters are as follows: 13.75 and 1100 mm focal lengths for OROCHI and TENGOO, respectively, based on an early mission plan document.These values differ slightly from the latest values of 3.23-13.52and 947.8 mm for OROCHI and TENGOO, respectively, as reported in Kameda et al. (2021).However, the differences do not affect our simulation results significantly.The pixel size of the charge coupled device (CCD) and the pixel number are the same in both cameras (5.5 µm and 3296 × 2472, respectively).For simplicity, the line-of-sight vector from each camera to the Phobos surface at each time point is assumed to be the same as that of LIDAR.The data interval is set to 300 s.In the actual mission, the imaging intervals by OROCHI and TENGOO during QSO operations might be around 3600 s due to on-board storage and data download capacity limitations.The effects of sparse sampling intervals on the estimation of the gravity field are described in Sect.4.5.
The coordinates on the image plane of each landmark captured by the on-board camera at each time are calculated as follows: to capture a landmark by the on-board camera, it must be 1) on the visible hemisphere as seen from the spacecraft, 2) illuminated by the Sun, and 3) the projection on the image must be within the camera field of view.We first compute two inner products, i.e., for a normal vector of a landmark and a vector from the spacecraft to the landmark, and for a normal vector of a landmark and a vector from the Sun to the landmark.These products are used for all landmarks defined on the Phobos surface for each observation time and we use them to select landmarks that satisfied Conditions 1) and 2).Subsequently, for the selected landmarks, projected coordinates on the image plane are calculated as follows: where X img , Y img is a projected landmark position onto the image plane, and r f and X pix , Y pix are the focal length and pixel size of the camera, respectively.The origin of the image plane is defined as the top-left corner, and (X 0 , Y 0 ) is an offset of the origin with respect to the image center.Moreover, λ ls and θ ls are the azimuth and elevation angles of the landmark direction, respectively, as seen from the spacecraft and calculated as follows: where x ls x ls , y ls , z ls is a landmark position in the space- craft-fixed frame and derived as follows: where x l is a landmark position in the Phobos-fixed frame, x is the spacecraft position in Phobos-centered inertial frame, and R 0 and R 1 are coordinate transforma- tion matrices from Phobos-fixed to inertial frame and from inertial to the spacecraft-fixed frame, respectively.Among the calculated landmark coordinates, we use those that satisfied Condition 3), i.e., 0 ≦ X img ≦ 3296 and 0 ≦ Y img ≦ 2972, as the observed landmark coordinates within the image plane at each observation time.Equations (6-8) allow the landmark coordinates on the image to be related to the spacecraft position, and thus, the data can be used to estimate the spacecraft orbit.
For random Gaussian observation errors with standard deviation 0.03 mm/s for 60-s integration time, 2 m and 0.5 pixel are added to all simulated observed Doppler values (range rate), LIDAR range, and landmark coordinates.In reality, LIDAR range observation error varies with altitude, showing values of 2 and 22 m at 100 m and 100 km, respectively (Senshu et al. 2021).However, since most of the error change with altitude is a biased component, we assumed an ideal case that the bias could be estimated and removed, and thus, we use 2 m as the standard deviation of the LIDAR range error in this study.As for the attitude control error of the spacecraft, which affects the line-of-sight direction of the LIDAR range and the image observations, random Gaussian errors with a standard deviation of 0.03° are associated with the lineof-sight vector at each time step.
For orbit and gravity field estimation, observation data are assumed to be available for seven consecutive Earth days without the orbital control for each QSO. (6) We assume that the Doppler observation data are available for 8 h a day acquired from any of the four stations, Usuda, Goldstone, Madrid, and Canberra.This is assumed to simplify the setting of Doppler observation times in the simulations.In an actual MMX mission, the Misasa station is mainly used for Doppler observations when the spacecraft is visible from Misasa.Observations at other stations are restricted during critical operations.Further, LIDAR range and image observations are assumed to be conducted when the spacecraft passes the dayside of Phobos when Doppler observations are not conducted.For the ascending trajectory, we assume that all Doppler, LIDAR range, and imaging observations are conducted from an altitude of 100 m above the Phobos surface to 50 km.

Simulation procedure
In the simulation, we assume homogeneous density of the Phobos model as the 'true' Phobos and estimate the Stokes coefficients of the gravity field using Doppler, LIDAR range, and landmark coordinate datasets acquired from QSO.We also assume two-layer internal structure density model as the 'true' Phobos, and estimate the inner layer density or the density boundary of the two layer using Doppler, LIDAR range, and/or landmark coordinate datasets acquired from the ascending trajectory.The initial state vector of each arc is also estimated for all simulations.Figure 3 shows an overview of the simulation flow followed in this study.
First, the initial state vector of the spacecraft and other estimated parameters, e.g., Stokes coefficients, are considered and their values are considered as the "true" values in the simulation.These values are used to generate a simulated, perfect spacecraft orbit by orbit integration and to calculate perfectly fitting observables.Second, to achieve a more realistic scenario, random observation errors are added to these observables leading to a set of "realistic" observables.Conversely, we add errors to the initial state vector of the spacecraft and other estimated parameters, and use them to generate a "computed" orbit and calculate "computed" observables from this orbit.Finally, we estimate the initial state vector and other estimated parameters that minimize the difference between "realistic" and "computed" observables by the leastsquares method.When dealing with multiple types of observables, weighted least-squares method is employed.We assume that each observation has independent measurement errors without cross-correlation and set the weight matrix , where σ i is the observation error for i th observation.As for the value of σ i , we use the standard deviation of the error for each type of observation, as described in Sec.2.3.The recovery precisions of the estimated solutions are evaluated by comparing them with the "true" values.In this study, we assume unbiased estimation for simplicity in the simulation.However, in the actual mission, various unpredictable factors, such as instrument failures or telecommunications troubles, may hinder ideal data acquisition for each observable, potentially introducing biases in the data.Consequently, there is a possibility that simple weighted least-squares methods may not yield unbiased estimates.Taking this into consideration, in the actual mission, it is necessary to conduct preliminary investigations by creating and testing numerous subsets of solutions to explore the contribution of different data types.
In this study, we assume that the shape model and the solar radiation pressure model are error-free, and thus, the parameters related to these models are not estimated.Moreover, we do not estimate the ephemeris assuming that the Phobos ephemeris is perfectly known.However, in Sect.4.4, we evaluate the effects of the error of the Phobos ephemeris on the Stokes coefficient estimations.In the actual mission, the Phobos ephemeris will be estimated and updated at the early stage of the Phobos observation phase using radio-tracking datasets, LIDAR, and a shape model.

Estimation of the gravity field coefficients from QSOs
Using the Doppler, LIDAR range, and landmark coordinate data observed from QSOs, we simulate the estimation of the initial state vector of each arc and Stokes coefficients of Phobos gravity field.In the actual MMX mission, spacecraft orbit control and/or attitude control is planned to be conducted at least once every few hours to a day.In consideration of this, the observation data are divided into 1-day segments and the initial state vector of each 1-day arc-length orbit is estimated as an arcspecific parameter.Subsequently, the Stokes coefficients are estimated as global parameters from the 7-day orbits.We select 7-day data to ensure that all each 1-day arc include not only Doppler observations, but also LIDAR and image observations, which are conducted exclusively when the spacecraft is illuminated by the sun.In the case of QSO-H, 3D-QSO-M, QSO-LA, and 3D-QSO-LA, the spacecraft orbits Phobos with orbital periods of approximately 32, 4, 1, and 1 day(s) in the inertial frame, respectively, half of which correspond to the duration of solar illumination.For QSO-H, we utilize consecutive 7-day data collected during periods of solar illumination for simulation.In the case of 3D-QSO-M, during consecutive 7-day periods, there are some 1-day arcs that are not illuminated by the sun, resulting in the unavailability of LIDAR and image landmark data.Therefore, for 3D-QSO-M, seven 1-day arcs containing Doppler, LIDAR, and image data are selected from the continuous 14-day observation data.In the case of QSO-LA and 3D-QSO-LA, as each 1-day arc inherently includes illuminated and nonilluminated periods, we employ data from arbitrarily selected consecutive 7-day arcs.A Phobos model with homogeneous density is assumed in this simulation.Random Gaussian errors having standard deviations of 10 m and 0.001 m/s are added to the initial spacecraft position and velocity of the true orbit for each arc, respectively.For the Phobos gravity field model, we assigned an initial error to each estimated Stokes coefficient.The error values are assumed to be on the order of 100%, while errors are not considered for the other nonestimated coefficients, while errors are not considered for the other non-estimated coefficients.Degree 1 terms are assumed to be defined values, and are not estimated.Figure 4, Table A1.1,A1.2, A1.3, and A1.4 show the estimated results of the gravity field.
As the orbit is 2D, estimating the Stokes coefficients other than the sectorial ones is difficult from QSO-H and QSO-LA.However, the other degree-2 term, C 20 , can also be estimated from QSO-H, because the orbital plane of QSO-H is slightly inclined with respect to the equator.The magnitude of the relative error depends on the orbital altitude and coefficient value.Table A2 shows the "true" Stokes coefficients and Fig. 5 shows their amplitudes.Comparing Figs. 4 and 5, we find that a coefficient with a large amplitude can be estimated better because its contribution to the total gravity acceleration is relatively large.From QSO-H (distance from the Phobos center: 100-200 km), C 20 , C 21 , C 22 , S 22 and S 33 terms are estimated with better accuracy than the 100% error given as the initial error.Particularly, the C 20 term is well-esti- mated with 0.0375% error.On the other hand, the S 21 term is poorly estimated compared to the other degree 2 terms.From the viewpoint of the orbital configuration of QSO-H, we do not observe that the S 21 term is par- ticularly difficult to estimate compared to other degree 2 terms.Rather, the difficulty in estimation is attributed to the fact that the value of S 21 is smaller than that of other degree 2 coefficients.C 21 , S 21 and S 22 terms are related to the orientation of Phobos' rotation axis.If the z-axis aligns with the polar principal axis of inertia, C 21 and S 21 become zero.On the other hand, S 22 becomes zero when the x-axis coincides with one of the equatorial inertia axes.In the case of Phobos, which is in synchronous rotation state around Mars, defining the x-axis towards Mars results in S 22 being zero.In this study, the origin and coordinate axes are directly utilized as defined in the shape model without adjustments to the origin or rotation axis directions.Consequently, the values of C 21 , S 21 , and S 22 are small; particularly, the S 21 term is smaller by an order of magnitude compared to the C 21 , and S 22 terms (Table A.2), but they are not entirely zero.From QSO-LA (distance from the Phobos center: 30-50 km), the sectorial coefficients up to degree 7, except S 77 term, are estimated with improved accuracy compared to the accuracy of the initial values.
From 3D-QSOs, in principle, not only sectorial but also full degree and order coefficients can be estimated owing to their latitudinal coverages.The estimation accuracy depends on the altitude and coefficient value.Although some small coefficients cannot be estimated with sufficient accuracy, estimating coefficients is possible up to approximately degree 4 from 3D-QSO-M and degree 5 from 3D-QSO-LA with improved accuracy than the initial coefficient, as shown in Fig. 4c, d. Figure 6 shows the degree variance of the difference between the model with Stickney's gravity anomaly and the homogeneous model, and the ones of the estimation errors for 3D-QSO-M and 3D-QSO-LA.From this figure, it can be observed that in the case of 3D-QSO-M, errors exceed the signal of the gravity anomaly at degrees 4 and above, whereas in the case of 3D-QSO-LA, even at degree 5, errors remain below the signal, allowing for the detection of the gravity anomaly.
The comparison of degree variances provides the advantage of understanding the resolution of detectable signals in the spectral domain, but it relies on globally developed spherical harmonic coefficients.In practical cases, even when it is expected to be challenging to detect spectrally resolved signals through comparisons of degree variances, it may still be possible to detect them if the amplitude of the local signal is sufficiently large.Therefore, we compare errors and signals in the spatial domain as well.Figure 7 depicts the estimation errors of gravity field in the spatial domain using up to degree and order 2, 3, 4 and 5.The errors for 3D-QSO-M are approximately ± 0.15, ± 2.65, ± 11.1, and ± 97.5 mGal for the truncated degree 2, 3, 4, and 5, respectively.Similarly, the errors for 3D-QSO-LA are approximately ± 0.02, ± 0.10, ± 0.39, ± 1.10 mGal for the truncated degree 2, 3, 4, and 5, respectively.
Figure 8 shows the global map of the regional density anomaly model described in Sect.2.1.These maps look similar to the one in Fig. 1, but it utilizes Stokes coefficients up to degree and order 2, 3, 4 and 5, respectively.The maximum amplitudes of the gravity anomaly of Stickney are 13.4,19.7,25.2 and 29.5 mGal,respectively. Comparing Figs. 7 and 8, in the case of 3D-QSO-M, the errors are smaller than the signal when truncated at degrees 2, 3 and 4. The signal-to-noise ratio is 2.3 when truncated at degree 4, which is slightly different from the expected result based on degree variance, as anticipated from the earlier description.When truncated at degree 5, the error is larger than the signal.On the other hand, utilizing data from 3D-QSO-LA allows for a clearer observation of Stickney's localized signal with higher resolution and a signal-to-noise ratio of 32, even when truncated at degree 5, enabling a more significant detection of the signal.Thus, considering the balance between the error magnitude and the spatial resolution of the signal, 3D-QSO-LA appears to be the preferable choice for detecting gravity anomalies at the scale of Stickney Crater than 3D-QSO-M, emphasizing the importance of observations with 3D-QSO at such a low altitude as LA.

Detection of differences in the internal density structure from the ascending trajectory
We first investigate whether distinguishing differences in the layered density structure models is possible by observations from the ascending trajectories.The spacecraft experiences a greater influence of gravity acceleration from mass elements in shallow regions closer to the surface because magnitude of gravity acceleration depends not only on the density, but also on the distance from each mass element.Therefore, the gravity acceleration on the spacecraft differs depending on the density structure in the radial direction.Figure 9a shows each trajectory of the MMX spacecraft assuming different Phobos density structure models, and Fig. 9b shows the differences in LIDAR range observations between the layered density structure model and the uniform density model for three different water contents.
The state vector of the MMX spacecraft is affected by the difference in gravity acceleration caused by the differences in the Phobos density structures.The acceleration differences are larger at lower altitudes and decrease with increasing altitude.However, the difference in position and velocity develops over time, reflecting the effect of the difference of the accelerations at low altitudes.After 10 min, the differences in the LIDAR range observations increase than the observation error of 2 m.Thus, detecting the differences in the layered density structure models is possible by LIDAR observations.The same simulations are also performed using landmark coordinate and Doppler observation data.For landmarks observed by the wide-angle camera OROCHI, even after 65 min, the coordinates differ by an average of 0.21, 0.47, and 0.80 pixels with respect to homogeneous models, for 10%, 20%, and 30% water contents, respectively.This is because gravity acceleration is dominant in the radial direction and is difficult to detect in image observations.For Doppler data, the sensitivity of the difference detection varies with time.High-sensitivity detection of the difference is possible when the position vector of the spacecraft as seen from the ground station is parallel to the velocity vector of the spacecraft.However, when the two vectors are nearly orthogonal, the sensitivity decreases.Thus, among the three observation types, LIDAR range observations are the most useful for detecting internal structures from the ascending orbits.
Next, we simulate the detection of regional density anomaly beneath the Stickney crater through observations from ascending trajectories.As shown in Fig. 10, the difference in LIDAR range values with respect to the homogeneous density model are larger than the observation errors after a certain period, indicating that the regional density anomalies can be detected by LIDAR range observation if the landing site is close to the anomalies.In this simulation, we detect positive regional gravity anomaly assuming that the impact causing Stickney Fig. 10 Differences in LIDAR ranges observed from the MMX ascending trajectory between using Phobos internal structure model with homogeneous density and with regional density anomaly beneath the Stickney crater.See Sect.2.1 for each Phobos internal density structure model Fig. 9 a MMX ascending trajectories assuming layered internal density structures of Phobos with 0, 10, 20, and 30% water content.View from the north pole of Phobos.b Differences in LIDAR ranges observed from the MMX ascending trajectory between using the Phobos model with layered internal density structure with 10, 20, and 30% water content, and with homogeneous structure (0% water content).See Table 2 for the Phobos internal density structure models resulted in a compression of the material lying below the crater (Richardson et al. 2002).However, detecting negative gravity anomaly with less dense area below the crater is also possible (Britt et al. 2002).In this case, the range difference shown in Fig. 10 would be negative.

Estimation of inner-layer density or density boundary of the layered density structure model from the ascending trajectory
We simulate the estimation of the density boundary radius or density of one of the two layers, which are parameters related to the layered density structure model, as well as the initial state vector from the ascending trajectory.Assuming Phobos has a two-layer internal density structure, in which its total mass M and total volume V are constants.As mentioned in Sect.2.1, we assume that the shape of the density boundary in the two-layer model is similar to that of the surface.In this case, common spherical harmonic coefficients can be used for both the inner and outer layers.We assume that the density boundary is located at a distance of i r times from the center to the surface, where 0 < i r ≦ 1 , and denote the densities of the inner and outer layers as ρ inner and ρ outer , respectively.Following from the law of mass conservation, The partial derivatives of Eq. ( 9) with respect to ρ inner and i r are as follows, respectively: and Thus, the partial derivatives of ẍ , which is the gravity acceleration of Phobos experienced by the spacecraft, with respect to ρ inner and i r are derived as follows: and Thus, parameters ρ inner or i r can be estimated through analysis of the spacecraft orbit.Because the two parameters ( ρ inner and i r ) degenerate, estimating both param- eters from dynamic orbit determination is impossible.To (9) M = ρ inner i r 3 V + ρ outer 1 − i r 3 V . ( estimate both parameters, other independent observation data, e.g., seismic tomography data, are required.In case of density estimation, we estimate ρ inner .However, estimating ρ outer as an estimation parameter instead of ρ inner is also possible.Once the density of one layer is esti- mated, that of the other layer can be uniquely determined by conserving mass and volume using Eq. ( 9).Moreover, in this simulation, we assume the 20% water content case model in Table 2 as the "true" Phobos density structure model for the simulations, and accordingly, the innerlayer density or density boundary is estimated.Similar to the estimation of Stokes coefficients from QSOs, random Gaussian errors having standard deviations of 10 m and 0.001 m/s are added to the initial spacecraft position and velocity.The initial parameter values with errors are set to 90% of the radius from the center for the density boundary radius and 1000 kg/m 3 for the inner-layer density, respectively.When only the Doppler dataset is used, the solution do not converge because of insufficient accuracy to detect the Phobos density structure, as mentioned in Sect.4.2.Moreover, when both Doppler and LIDAR range data are used, the parameters cannot be estimated well.This is because the LIDAR footprint position on the Phobos surface cannot be specified from the line-of sight laser range data only, and correcting the wrong footprint position to the right one is difficult if the "computed" trajectory differs considerably from the "true" trajectory (Yamamoto et al. 2020).Conversely, when the Doppler and landmark coordinate data acquired from OROCHI are used, solutions can be obtained.The relative errors of the estimated results with respect to the "true" values are 3.15% and 1.20% for the radius of density anomaly and innerlayer density, respectively.Further, better results can be derived when LIDAR range data are used along with the Doppler and landmark data, because LIDAR range observations are sensitive to the radial-direction changes in the gravity field caused by the layered density structure of Phobos.In this case, the relative errors are 0.220% and 0.0541% for the radius of density anomaly and inner-layer density, respectively.

Effects of Phobos ephemeris errors
In the simulations described in Sects.4.1-4.3,we assume that no error is contained in the Phobos ephemeris.However, estimating the effect of ephemeris errors on parameter estimation is a concerning issue.This section describes the simulation of the Stokes coefficient estimation using the 3D-QSO-LA data when the Phobos ephemeris is inaccurate.
We assume an error existed in the time of perigee passage of the Phobos orbital elements, which is one of the orbital elements that can be determined only inaccurately by astronomical observations.An error of 0.01 s is considered for the perigee passage of the osculating orbital element at the initial epoch.This error value corresponds to the order of 10 m and 0.001 m/s for the initial position and velocity, respectively, in the Cartesian coordinates.Further, we estimate Stokes coefficients of the Phobos gravity field up to the degree and order 5, including degree-1 coefficients, which correspond to the center of mass coordinates and are likely sensitive to the ephemeris error.
Figure 11 and Table A.3 show the simulation results.The accuracy of the estimated coefficients deteriorates compared to that shown in Fig. 4d, where degree-1 terms are defined and not estimated.Moreover, the error of the perigee passage epoch affects not only the estimation accuracy of degree-1 coefficients, but also of higher coefficients.Therefore, it is important to use an accurate ephemeris, and in the actual mission, simultaneous estimation of it is crucial.This is expected to be implemented in the early stages of Phobos observation, as described in Sect.3.

Effects of sparse sampling interval of the image landmark coordinate data
Figure 4 presents the results of gravity field estimation from QSOs using LIDAR range, Doppler, and image landmark coordinate observation data.The sampling interval is set to 300 s for image landmark data.However, the actual acquisition interval during the QSOs might be approximately 3600 s to save the Fig. 11 Relative errors of the estimated Stokes coefficients of the Phobos gravity field with respect to the true values, derived using observation data from 3D-QSO-LA.Results for the case when Phobos ephemeris has an error of 0.01 s in the time of perigee passage Fig. 12 Relative errors of the estimated Stokes coefficients of the Phobos gravity field with respect to the true values, derived using observation data from 3D-QSO-LA.Results for the case when the sampling interval of image landmark data is 3600 s.Other conditions are the same as in Fig. 4d data-downloading capacity, as mentioned in Sect.2.2. Figure 12 and Table A4 show the estimation results of the gravity field from 3D-QSO-LA when the sampling interval of the image landmark data is 3600 s; even with a sparse sampling interval, the accuracy is comparable to the results represented in Fig. 4d.

Discussion
Based on the results given in Sect.4, we discuss how information of the internal density structure of Phobos can be improved by the MMX mission at global, regional, and local scales.As stated by Le Maistre et al. ( 2019), a large-scale internal density structure can be evaluated from MOI.This approach is especially useful for radially layered density models that are simpler than the method used to evaluate complex pattern differences of the Stokes coefficients up to higher degree and order, which was also proposed in Le Maistre et al. (2019).To calculate MOI, libration amplitude of Phobos is required in addition to the degree-2 Stokes coefficients, C 20 and C 22 .Matsu- moto et al. (2021) assumed that Phobos has a two-layered internal density structure of the density boundary which is 0.9 times of the Phobos radius, and that 10% of water ice mass fraction is localized in the outer layer.They showed that to detect this layered density structure, MOI with 3% accuracy is required, and thus, the libration amplitude and degree-2 Stokes coefficients should be determined with 2%-3% accuracy by MMX observations.In the MMX mission, the libration amplitude accuracy is expected to improve significantly by the high-resolution optical observations.As shown in Fig. 4, for the Stokes coefficients, the required accuracy is sufficiently achievable for C 20 and C 22 by QSO-H observations.Further- more, QSO-LA observations improve the accuracy of C 22 , although they could not determine C 20 .Thus, QSO-H and QSO-LA observations together can be used to derive the accurate value of MOI, from which largescale radial density structure can be constrained.
The center of figure and the center of mass of Phobos are obtained from the shape model and degree-1 Stokes coefficients, respectively.Both values are expected to be greatly improved in the MMX mission and the difference in the two values can be used to detect the dichotomic distribution of the internal mass.The results described in Sect.4.4 show that the accuracy of the Phobos ephemeris affects the estimation of degree 1 terms.This is because the ephemeris represents the motion of Phobos' center of mass in the inertial frame.Therefore, estimating Phobos ephemeris simultaneously in the actual mission is important.
Using the datasets from 3D-QSO observations, we can estimate not only the sectorial terms, but also the full-matrix gravity field coefficients, which are useful for estimating the internal density structure, especially for identifying regional density anomaly.Le Maistre et al. (2019) introduced an overlapping factor to evaluate whether distinguishing different internal density structure models of Phobos from each other is possible from the corresponding Stokes coefficients, and reported that the homogeneous density model, two regional density anomaly models below the Stickney crater, and two-layered density models clearly differed in the coefficients at the specific order and degree, which indicates that these models can be significantly distinguished.Although they used Stokes coefficients up to degree and order 10 for this evaluation, their results showed that the difference can be distinguished even up to degree and order ~ 5, which can be detectable by 3D-QSO-LA observations of the MMX mission, as described in Sect. 4. Conversely, especially for regional density anomalies, the anomaly location is more directly and easily identified using the gravity anomaly map rather than the spectral coefficients (Fig. 8).In the spatial domain, geographical correlations of density anomalies with topography and other remote sensing data can be found more directly, which may help understand the origins of the density anomalies.Information on the magnitudes and locations of density anomalies is also useful for discussing scientific interest and technical landing safety when selecting landing sites.In the MMX mission, the observation from QSO-LC (± 20 × ± 27 km), which is lower QSO than (3D-) QSO-LA, is also planned.Higher-degree and order Stokes coefficients are expected to be determined using observational data from the lowaltitude QSO.However, in such a low-altitude orbit, the spacecraft trajectory is more sensitive to the gravity field of Phobos.Hence, if the initial gravity field model has a large error, the difference between the computed and true orbit increases, making estimations of orbit and gravity field coefficients difficult.Therefore, the lowdegree coefficients estimated at high altitudes should be updated from the initial coefficients before estimating the orbit and gravity coefficient using low-altitude data.
Various types of observation data acquired from QSOs at high, middle, and low altitudes are yet to be comprehensively analyzed and used to select landing sites on the Phobos surface for sample acquisition.By then, we will be able to provide information on the magnitude and location of regional density anomalies with spatial resolutions of several kilometers in half-wavelength on the Phobos surface, besides the large-scale radial density structure obtained from MOI.Although locations with significant density anomalies are not necessarily selected as landing sites of the MMX spacecraft, the difference between the actual and predicted ascending trajectory after landing can help determine find local density anomalies around the landing site.This information may contribute to the study of the origin of the acquired samples.As shown in Sect.4.3, estimating the inner (or outer) layer density or density boundary is possible using the observation data from the ascending trajectory, assuming a layered internal density structure of Phobos.This is basically equivalent to the evaluation of radial density structure by MOI, but if both results are not consistent, it might reflect a local difference in the density or boundary depth.If significant density anomaly is expected near the landing site from the global gravity model created by QSO observations, we can estimate either the density or the anomaly depth by assuming a 2D location of the density anomaly.
As mentioned previously, the use of accurate Phobos ephemeris is extremely important in these gravity field estimation processes.Additionally, using an accurate shape model for orbit and gravity field estimation is important.Moreover, accurate orbit and gravity field models contribute to the improvement in the shape model in terms of providing accurate camera image positions.As stated in Matsumoto et al. (2021), these aspects will be improved by iterating the orbit and gravity field estimation and shape model determination during the mission.

Conclusion
The degree-2 Stokes coefficients of the Phobos gravity field, estimated from the QSO-H and QSO-LA observation data will be used for the calculation of MOI, along with the libration amplitude, which is also estimated in the mission by optical observations.The accuracy of the estimated degree-2 coefficients is sufficient to meet the MMX mission requirement related to MOI.Moreover, 3D-QSO-LA is not nominal in the current mission plan, but if it is realized, the Stokes coefficients up to degree and order 5 could be estimated from the observations at this orbit.This spatial resolution and accuracy are sufficient to detect a large density anomaly beneath the Stickney crater expected in one of the large-scale internal structure models of Phobos.The large-scale internal density structure of Phobos can be inferred by the orbit and gravity field determinations from QSOs at different altitudes and can be used to constrain the origin of the body.Additionally, the detection of gravity anomalies from ascending trajectories may help in interpreting the origin of the return samples.Thus, the gravity field derived in the MMX mission is expected to provide better insights into the density structure of Phobos at global, regional, and local scales.
In this study, we conducted simulations mainly focusing on the estimation of gravity field and internal density structure of Phobos.However, in the actual MMX mission, it is also important to estimate other geodetic parameters, e.g., ephemeris, rotational parameters, the tidal Love number k 2 , and so on.This would contribute to a comprehensive understanding of the internal structure and evolutionary process of Phobos.

Fig. 1
Fig.1Gravity anomaly map of the Phobos internal structure model with regional density anomaly beneath the Stickney crater with respect to that with homogeneous density depicted using Stokes coefficients up to degree and order 180 (see Sect. 2.1 for model development)

Fig. 3
Fig. 3 Schematic diagram of the simulation flow

Fig. 4
Fig. 4 Relative errors of the estimated Stokes coefficients of the Phobos gravity field with respect to the true values, derived using observation data from a QSO-H, b QSO-LA, c 3D-QSO-M, and d 3D-QSO-LA

Fig. 6
Fig. 6 Degree variance of the difference between the model with Stickney's gravity anomaly and the homogeneous model, and the ones of the estimation errors for 3D-QSO-M and 3D-QSO-LA

Fig. 8
Fig. 8 Same as Fig. 1, but using Stokes coefficients up to degree and order a 2, b 3, c 4 and d 5

Table 2
Homogeneous and two-layer Phobos internal density structure models used for the simulations described in Sect.4.2