Covariance analysis of periodic and quasi-periodic orbits around Phobos with applications to the Martian Moons eXploration mission

Understanding the internal composition of a celestial body is fundamental for formulating theories regarding its origin. Deep knowledge of the distribution of mass under the body’s crust can be achieved by analyzing its moments of inertia and gravity field. In this regard, the two moons of the Martian system have not yet been closely studied and continue to pose questions regarding their origin to the space community; thus, they deserve further characterization. The Martian Moons eXploration mission will be the first of its kind to sample and study Phobos over a prolonged period. This study aims to demonstrate that the adoption of periodic and quasi-periodic retrograde trajectories would be beneficial for the scientific value of the mission. Here, a covariance analysis was implemented to compare the estimation of high-order gravitational field coefficients from different orbital geometries and for different sets of processed observables. It was shown that the adoption of low-altitude non-planar quasi-satellite orbits would help to refine the knowledge of the moon’s libration angle and gravitational field.


Introduction
The origin of the Martian moons, Phobos and Deimos, remains unknown to the astronomical community. It is not known whether they are asteroids from the Main Belt trapped by Mars' gravitational field or the result of a giant collision between proto-planets [1,2]. However, through analysis of the moons' composition, hints pertaining to their origin could be revealed. In this regard, studying the moons' libration motion and gravitational field could provide insight into the distribution of mass under their crusts [3,4].
The Martian Moons eXploration (MMX) mission, a JAXA project scheduled for launch in the mid-2020s, could play a crucial role in exploring the two Martian moons. It will be the first mission to target these moons since the Phobos 1 and 2 missions, which were unfortunately lost before reaching the mission targets [5,6]. After an initial insertion into the orbit around Mars, the spacecraft will be moved to a set of quasi-satellite orbits (QSO) around Phobos for a series of proximity operations, including surface landings to gather samples.
The QSO orbiting phase will provide valuable data for analyzing Phobos's dynamics and gravitational field [7]. At the time of writing, planar QSOs with different reference altitudes around Phobos are under consideration for further characterization and study of the geophysical environment of the Martian moon. In particular, five planar QSOs were selected for the observation of the satellite, as reported previously [8].
However, considerations and preliminary navigation analyses emerging from the Centre national d'études spatiales (CNES) revealed that the selection of lowaltitude quasi-periodic orbits could be beneficial for further characterizing Phobos's dynamical environment. As part of these preliminary investigations, mid-altitude spatial retrograde orbits (3DQSOs) and swing QSOs have been proposed to enable global coverage and refine the gravitational field of the small body [8,9]. This study elaborates on the original findings of JAXA and CNES to further investigate the navigation requirements for lowaltitude periodic and quasi-periodic trajectories around Phobos. In-plane and out-of-plane deviations from the proposed MMX QSOs are considered to achieve a better estimation of the moon's state and gravitational field coefficients.
A covariance analysis tool enabled by unscented Kalman filters (UKF) was implemented to analyze the spacecraft's motion along periodic and quasi-periodic orbits around Phobos. A mid-fidelity model that considers the roto-translational effects of an oblate Mars and a second-order, second-degree gravitational field of Phobos has been implemented; moreover, in this new model, the previously mentioned reference trajectories have been designed and optimized. The covariance analysis works by processing the measurement data types that will be available throughout the mission [10]. Practically, MMX's navigation is designed to rely on ground-based observables, that is, range and range rate measurements from the Deep Space Network (DSN) antennas, together with Lidar measurements and line of sight (LOS) direction vectors from the onboard cameras.
The remainder of this paper is divided as follows. Section 2 describes the adopted dynamic model, which includes the translational and rotational dynamics of an ellipsoidal moon around Mars. Section 3 outlines the procedure for extracting reference orbits in the new dynamic model. Section 4 describes the observables' models, followed by a description of the tool designed to perform the covariance analysis in Sections 5 and 6. Section 7 presents the results of the analysis, focusing on the reconstruction of Phobos's libration angle, moments of inertia, and gravitational field.

Dynamic model
Consider the motion of MMX, as seen from the Mars Mean Equator and IAU vector of the J2000 frame (MARSIAU), as shown in Fig. 1. The state vector of the spacecraft is hereby referred to as X MMX = [r, v], where r = [x, y, z] and v = [ẋ,ẏ,ż] are the position and inertial velocity vector of the spacecraft as seen from Mars, respectively. The state of Phobos is modelled in polar coordinates, assuming that motion occurs on the equatorial plane of Mars. As previously described [11], the coupled dynamics of the Mars-Phobos system can be captured by an 8-dimensional state that includes the radial distance of Phobos from Mars r Ph , Phobos's true anomaly θ, Mars and Phobos's libration angles (Φ M and Fig. 1 Schematic representation of the problem geometry with MMX's and Phobos's state vector components. Some of the reference frames involved in the analysis are also represented (the frame rotating with the moon's true anomaly, the Phobos-fixed reference frame, and how are they defined with respect to the perifocal frame. Φ Ph ), and their respective time derivatives. In addition, to analyze Phobos's internal structure, its moments of inertia are included in the state (I Phx , I Phy , and I Phz ). Lastly, some biases affecting the observables will be monitored (ϵ ρ , ϵρ, and ϵ Lidar ), yielding the 20-dimensional state vector of Eq. (1): The time evolution of the problem's state vector X is described by the set of nonlinear first-order differential equations in Eq. (2): where a M denotes the acceleration produced on the spacecraft by Mars, a Ph is that caused by Phobos, and a PhM is the gravitational transport between Mars and Phobos affecting MMX. The accelerations acting on the spacecraft and the time derivative of the moon's state vector components are described separately in the following subsections. Moreover, Phobos's moments of inertia and the observables' biases are constant during the analysis; therefore, their time derivative is null.

Phobos dynamics
The moon's motion is described with a planar full 2-body problem considering the coupling between the rotational and translational motion of the two bodies [11]. It is described with respect to the moon's orbit perifocal reference frame using eight ordinary differential equations, which are derived from the mutual potential V : As shown in Eq. (3), the system's motion is influenced by the reciprocal gravitational attraction and the rotational dynamics of the two main bodies. They are modelled as an oblate main body (Mars) and an ellipsoidal secondary body (Phobos) [11]. The equations of motion for the system are where ν is the main body mass fraction, and the partials of the potential with respect to the radial distance and Phobos libration angle are These dynamics can be propagated in the Mars-Phobos case and compared to the corresponding data extracted from NASA's SPICE Toolkit [12]. This comparison was performed using a two-day time window, and as shown in Figs. 2 and 3, the model has been shown to properly describe the behavior of the planet-moon system over time.

MMX dynamics
The dynamics of MMX are modelled as an N -body problem that includes the gravitational influence of Mars and Phobos on the spacecraft. Mars's influence on MMX is composed of the Keplerian acceleration and its oblateness gravitational terms. Therefore, it was computed directly in the MARSIAU reference frame as where µ M is the red planet's gravitational constant, R M is Mars's mean radius, and r is MMX's position vector magnitude.
As previously stated, Phobos is modelled as a triaxial ellipsoid whose gravitational influence can be modelled according to the gravitational potential described in Ref. [13]: where δ is the latitude, λ is the longitude of the spacecraft with respect to Phobos's principal body-fixed reference frame, Phobos's reference sphere has radius r 0 Ph , and the distance of the spacecraft from its center is |r − r Ph |. To recover this vector and calculate a Ph , three rotations are required for MMX's position vector from the MARSIAU reference frame to Phobos's principal axis frame. The resulting acceleration vector must then be rotated back to the MARSIAU reference frame for propagation.
To complete the formulation of the N -body problem affecting MMX's motion, the transport term consisting of the attraction to Mars caused by Phobos (a PhM ) must be derived. Following the third principle of dynamics, a PhM corresponds to the acceleration caused by Mars on Phobos once scaled by the mass ratio of the two celestial bodies.
The acceleration impressed by Mars on Phobos can be derived from the mutual gravitational potential [11]. The only complication might lie in deriving the acceleration components in the MARSIAU reference frame because the potential reported in Eq. (3) is a function of eight different states. Nevertheless, Phobos's radial distance from Mars and its true anomaly can be viewed as spherical coordinates once in the perifocal reference frame. The partial derivatives of the potential with respect to these parameters can be obtained, and the partial derivative with respect to the radial distance is expressed in Eq. (5) and the one with respect to the true anomaly is reported in Eq. (9) together with the partial with respect to the polar angle φ, which is null, being Phobos's motion planar in the perifocal reference frame.
As a result of Eqs. (8), (9), and (10), the acceleration impressed by Mars on Phobos can be recovered by exploiting the differentiation chain rule. The partial derivatives of the spherical coordinates with respect to x Ph , y Ph , and z Ph are known and can easily be derived. The acceleration described in Eq. (8) can be derived and rotated to the MARSIAU reference frame to be included in MMX's motion.
The resulting dynamic model is more complex than that adopted previously in Ref. [8], in which the mission reference orbits were identified. The identified trajectories were defined in a Hill problem centered on a tidally locked Phobos, meaning that the libration movement of the moon was completely neglected. Moreover, the influence of Mars's oblateness on the spacecraft's motion was not considered. Therefore, in Section 3, a transition into these trajectories in our new dynamic model is investigated.

Reference trajectory definition
The existence of the periodic and quasi-periodic orbits, which are the targets of the analysis, is not guaranteed in the dynamic model set up for this analysis. Baresi and Kawakatsu [8] originally generated distant retrograde orbits, also known as quasi-satellite orbits, for the Mars-Phobos system by using predictor-corrector techniques in the circular Hill problem (CHP) for a tidally locked Phobos. Five relative trajectories were selected: a highaltitude orbit (  The same predictor-corrector techniques are not straightforwardly applied in the new dynamic model that includes the libration angle of Phobos, and an alternative method to design trajectories in this system is required. In Ref. [8] the authors transitioned the preliminary trajectories from the CHP to the fullephemeris model of the Mars-Phobos system. A similar procedure was followed in this analysis to identify reliable initial conditions for the propagation of the analysis's target trajectories, as described below.
Periodic planar QSOs were first calculated in the CHP and later re-propagated for 20 days to store the periapsis and apoapsis points, where These N nodes are characterized by different time stamps X MMXi (t i ) separated by half of the orbit's period. To identify the same periodic trajectories that move following the new dynamics, a few more passages need to be taken: First, Phobos is propagated following the new dynamics [11] and its states at the corresponding time of the nodes are stored. The state vectors of the nodes of MMX's Hill problem are translated to the Phobos-fixed frame. These initial conditions are then rotated to the MARSIAU reference frame to be used as initial estimates for a continuation procedure that consists of patching forward and backward branches propagated using the newly defined dynamic MMX model, as shown in Fig. 5.

Phobos Phobos
Backward propagation F o r w a r d p r o p a g a t i o n Defining the MMX state vector at the i-th node as X MMXi = [r i , v i ], the cost function is minimized as Eq. (13): Nonlinear constraints are added to force the patching points to actually patch, together with constraints related to the continuity of time between the branch propagations. Moreover, constraints on the radial velocity of the nodes are added to facilitate convergence, and the initial and final nodes are bound not to move outside of a box of dimensions [k x , k y , k z ]: The trajectory optimization procedure then adjusts the nodes and epochs until the trajectory's continuity is guaranteed. If the procedure determines a candidate orbit, the cost function defined in Eq. (13) converges to zero. Practically, a tolerance value for the cost function must be set and considered to be acceptable as zero. An overall ∆V Tot ⩽ 10 −6 km/s was considered satisfactory. An example of the newly found QSOs can be found in Fig. 7.
A similar procedure was applied to transition the 3D and swing QSOs at five different orbit heights. The main difference is the definition of the initial guess for the continuation procedure, which is similar to that described in detail for the transition into a full-ephemeris model [8]. Both the 3D and swing QSOs can be calculated in the CHP using numerical procedures, as previously described [14]. Both of these orbits' families densely cover  the surface of 2D invariant tori. 2D quasi-periodic tori are characterized by two incommensurate frequencies (ω 1 and ω 2 ), as described previously [8]. Torus maps can be generated to map trajectories in R 6 onto trajectories in torus coordinates and vice versa. Figure 6 shows an example of a torus map that reports the radial velocity as a function of the torus coordinates for the 3DQSO-La orbit. Interestingly, for quasi-periodic trojectories (QPTs), unlike in the periodic orbit case, these points are not always aligned with the Mars-Phobos direction. As described previously [8], it is possible to extract the zero-radial velocity contours (red lines in Fig. 6) and calculate the spline coefficients to detect the intersections between the red curves and arbitrary quasi-periodic trajectories. The latter passage is the key to creating a time series of nodes along the surface of the considered quasi-periodic invariant manifolds to be used in the transition procedure to the new dynamic model. Once the nodes are available, they can be translated into the new dynamic model using the same optimization procedure outlined for planar QSOs. In the case of swing QSOs, the nodes are bound to remain in the equatorial plane and do not deviate from their original CHP position with respect to the moon by more than a few kilometers. Figures 8 and 9 show the final output of this continuation procedure, validating the strategy to generate 3D and swing QSOs around Phobos in the new model. Figures 7,8,and 9 show that the moon's irregular gravitational field, coupled with its libration motion, has  a disrupting effect on the periodicity of the orbits, which is particularly important in the vicinity of Phobos.
Once the reference trajectories are available, they can be used to finally create the observables for the "real" spacecraft to be processed in further analyses. The mathematical models used to describe them are reported in Section 4.

Measurement model
Four orbit data types were considered for this numerical analysis and are processed by the navigation system [10]: range and range rate from Earth-based ground stations, Lidar measurements spanning the moon's surface, and LOS direction vectors towards the natural satellite. The mathematical model and acquisition frequency for each are described in the following subsections.

DSN-based observables: Range & range rate
The first two observables simulated and processed in these analyses are the range and range rate measurements between Earth-based ground stations and MMX. Range measurements are typically obtained by measuring the time elapsed between the emission and reception of a signal from an antenna of the Deep Space Network to the spacecraft, or vice versa. Range rate measurements are extracted from a measure of the Doppler effect caused by the relative radial velocity between the emitter and receiver antennas. The covariance analyses carried out in this study used an idealized range and range rate measurement between MMX and the Earth, thereby delegating clock errors, atmospheric effects, and time delays to the bias ϵ ρ [15]. Following these assumptions, the range measurements are defined as where r E and r Sti are the position vectors of MMX and the ground station, respectively, as seen from the Earth-centered J200 reference frame, and ϵ ρ is the range bias, which is also considered in the covariance analysis.
whereṙ E andṙ Sti are the velocity vectors of the spacecraft and ground station, respectively, with respect to the Earth-centered J200 reference frame. Note that occultation checks were performed throughout these numerical simulations, thereby preventing measurements from being collected whenever MMX was below the stations' local horizons or eclipsed by Mars and Phobos. The DSN antennas in Madrid, Canberra, and Goldstone have been considered, and measurements are collected at intervals of once per hour for the range and once per ten minutes for the range rate [16].

Lidar
Lidar is critical for planetary and small-body missions. A Lidar device typically emits an electromagnetic signal and listens for the echo of the same signal once it bounces off of the surface of a celestial body. As an electromagnetic signal, it travels at the speed of light; therefore, measuring the time distance between the emission and the reception of the echo provides information on the relative distance between Phobos and MMX. Considering the ellipsoidal shape of Phobos as a first-order approximation yields a state-observation relationship of the form [17]: where R Ph is the radius of the ellipsoidal Phobos as a function of the spacecraft's latitude (λ) and longitude (ϕ), defined with respect to the moon-fixed reference frame. As per Ref. [18], MMX's Lidar has an acquisition frequency ranging between 1 and 4 Hz. An analysis has been performed regarding the effects of this frequency on estimation performance, and it was found that there is a plateau for more than one measurement every 2 min. Therefore, this was identified as the Lidar acquisition frequency for the analysis. MMX's onboard Lidar has an operational range of 100 m-100 km; therefore, no measurements outside this range are collected.

Line-of-sight (LOS)
The MMX mission will be equipped with two navigation cameras for exploration purposes which will be used in the proximity phase to characterize the red planet's moons [19]. The TENGOO camera is a narrow-angle camera that can be used to identify Phobos's geomorphological features by pointing towards the spacecraft's nadir [20]. Owing to this pointing requirement, the camera can be used to generate LOS measurements from the spacecraft to the center of the moon.
Let us define a rotating coordinate system centered on the spacecraft and oriented such that the first axis, x CAM , is perpendicular to the focal plane of the camera and anti-parallel to the relative position vector of MMX as seen from the center of Phobos;ẑ CAM is parallel to the cross product betweenx CAM and the third axis of the J2000 frame, namelyẐ;ŷ CAM completes the righthanded triad: If the unit axes of the camera frame are provided, it is possible to calculate the projections of the relative position vector along the second and third components of the camera frame [21]: where η and ζ are the components of the relative position vector along the second and third axes of the camera frame, respectively. These are the observables for covariance analysis. A picture is collected and processed every hour for navigation purposes unless Phobos is in an eclipse, and therefore dark. Pictures are also discarded whenever the spacecraft is in opposition to the Sun, as observed from the moon: In other words, if the dot product between the vector from Phobos to MMX and the sunlight direction is positive. With all observables available, all that remains to be defined for the analysis setup is a tool for processing and retrieving the problem's state vector uncertainties. The setup of this covariance analysis tool is defined below.

UKF-based covariance analysis
A covariance analysis tool based on unscented Kalman filter theory [22] has been derived. 2N sigma points (where N is the state dimension) are commonly chosen when UKFs are adopted [22,23], based on a square-root decomposition of the a priori covariance matrix. These sigma points propagate through the nonlinear dynamics given by Eq. (2), and then the covariance is derived from their final distribution.
A few quantities must be defined to provide a detailed description of the tool. The first is {W i }, a set of 2N scalar weights defined as [22]: where λ = N (α 2 − 1) and α is a scaling parameter that guides the spread of the points around the mean of the distribution and is usually included in the interval 10 −4 ⩽ α ⩽ 1. For this analysis, α was set as close to 1 as possible in order to represent the original distribution of the sigma points around the reference trajectory as effectively as possible, ensuring that each is propagated correctly and does not follow Phobos. Second, the a priori state and covariance matrix are required. Here, we define the notation E[X], which represents the expected value of an arbitrary quantity X, andX, which represents the mean value of that quantity's distribution; the a priori state and covariance matrix of the state vector at t 0 are The covariance analysis tool propagates the reference trajectory, starting fromX 0 , through the dynamics for the duration of the analysis. For any covariance analysis tool, the estimated state is set to match the reference trajectory generated by the trajectory optimization procedure described in Section 3. Given the a priori covariance matrix P 0 , N initial sigma points can be generated as where η = √ N + λ and propagates until the next observation epoch. Generalizing the k-th step of the analysis, the time update of the covariance is performed at every epoch as where Q v denotes the process noise covariance. Only the stochastic acceleration acting on the spacecraft is considered; therefore, it is reflected in the described dynamics in terms of process noise covariance, as [15]: with Q = σ 2 u × I (3×3) and Γ = [∆t 2 /2 × I (3×3) ; ∆t × I (3×3) , 0 (14×3) ] where ∆t is the propagation time between observations. This results in a process noise covariance matrix with terms on the diagonal up to the 6th element, and extra diagonal terms only in the spacecraft positions and velocity coupling.
The covariance's time update is used to derive a new set of sigma points at time t k , which are used to derive the innovation covariance P yy and cross-covariance matrix P xy : where G[X k|k−1 ] represents the state-observation relationship (i.e., any of Eq. (15), (16), (17), or (19) depending on the measurement type under consideration),ȳ k is the observation vector derived from the reference trajectory, and R n is the measurement noise covariance. Using Eq. (24) the optimal Kalman gain and the update of the covariance matrix at the k-th step can be computed as P k is the estimated covariance at the k-th step and is used in the (k+1)-th step to generate the sigma points that are propagated along the reference trajectory. This process is repeated for all observations available throughout the duration of the covariance analysis.

Analysis setup
The main drivers of this analysis in the first instance were the covariance envelopes of the MMX state vector. The root mean square (RMS) values of the MMX's position and velocity vector component uncertainties were used as figures of merit to quantify the estimation chain performance.
where σ x , σ y , and σ z are the position vector component uncertainties of the spacecraft, and σ vx , σ vy , and σ vz are the velocity vector component uncertainties.
The second parameter of interest was the libration angle of Phobos. If properly observed, the moon's libration motion could be key to recovering information about its internal composition, as it is directly related to its moments of inertia and the gravitational torque exerted by Mars on Phobos [10,24].
The behavior of the uncertainties related to the moon's moments of inertia (MOI) is also analyzed because of their direct link with its gravitational field. In particular, the gravitational harmonic coefficients C 20 and C 22 can be recovered as functions of the MOI, as [25]: The uncertainties in these gravitational harmonic coefficients can be derived directly from the MOI uncertainties as [26]: where σ 2 Iii is the standard deviation of the MOI with respect to the i axis, and σ IiiIjj is the covariance between I ii and the MOI with respect to the j axis. None of these quantities have been accurately determined for either of the Martian moons in the state-of-the-art. Nevertheless, geodetic measurements of the gravity field and the rotation parameters of a body provide constraints on the moons' internal structure, reflecting their origin and evolution, which is one of the main targets of the Japanese mission [10].
Observations between MMX and the DSN ground stations are generated together with Lidar measurements of the MMX-Phobos distance and LOS measurements, while the spacecraft moves along the orbits and is updated with the different frequencies that were previously stated. All of the measurements are corrupted with zero-mean Gaussian white noise, whose standard deviations are reported in Table 1. The white noise values for the ground station data were obtained from Ref. [16], and the noise added to the Lidar measurements was recovered from technical data on MMX's onboard instrumentation [18]. In the same way, the white noise added to the camera's LOS direction vector was derived from TENGOO's specifications, and 1σ was considered equal to the camera's field of view divided by the number of pixels in the sensor. A stochastic acceleration of 10 −10 km/s was used as the process noise in the covariance analysis tool to consider the solar radiation pressure (SRP) and non-gravitational accelerations (NGA) affecting MMX's motion.
The a priori uncertainties for the states at the beginning of the covariance analysis are listed in Table 2. MMX's initial uncertainties were assumed, as performed previously [7], to have 3σ ≈ 1 km for the position vector components and 3σ ≈ 1 m/s for the velocity vector components. Large initial uncertainties in Phobos's libration angle were selected to demonstrate the observation capacity of the estimation chain. Those related to the moon's MOI were recovered from the C 20 and C 22 uncertainty values [10] using the uncertainty propagation rule described in Section 5.
One-week data acquisition campaigns were simulated for each orbit, and the data were fed into the covariance analysis tool to study the state vector components' 3σ covariance envelopes' behavior.

Comparison between orbit heights
It is interesting to observe the quality of the information that can be achieved by moving between different planar periodic QSOs and the altitude from the moon at which different features of the orbital environment can be observed.
As shown in Fig. 10, there is no significant difference in the position and velocity vector uncertainties that can be reconstructed from the covariance analysis of the five baseline QSO trajectories. The RMS of the 3σ envelopes for the position vector components ranges in 10 m ⩽ 3σ x ⩽ 50 m and the velocity vector components' RMS is in the range of a few millimeters per second, with QSO-Lc being the least accurate.
Regarding the moon, the same conclusions as those for  MMX can be stated: Phobos's state vector components can be characterized with the same degree of accuracy for all five orbits. The radial distance from Mars has a 3σ on the order of a few millimeters, and its true anomaly uncertainty falls below 10 −4 deg. The libration angle shows uniform behavior between the different heights that are directly observable from the Lidar. Its uncertainty oscillates in the range of 5 × 10 −2 deg ⩽ 3σ Φ Ph ⩽ 10 −1 deg, which provides valuable information for the reconstruction of the body's internal composition [10].
In general, the more closely the spacecraft approaches the moon, the better the final estimation of the gravitational field; however, there appears to be a limit to the achievable gravitational coefficient uncertainties, which depend on the distance at which the spacecraft orbits the moon and the amount of process noise adopted in the simulation. The larger the distance at which MMX orbits Phobos is, the weaker the moon's irregular gravitational field is, and the process noise hides part of this dynamic. This creates a bottleneck value for the observability of the C 20 and C 22 coefficient uncertainties that follows from the lack of observability of the moments of inertia. In particular, there seems to be a limit to the capability of the spacecraft's periodic motion to properly observe the two moments of inertia on the moon's equatorial plane, I xx and I yy . This aspect and the capability to better observe the I zz moment of inertia, an area with potential for improvement, are discussed in Section 7.2, where the effect of the orbit geometry is analyzed.

Comparison among processed measurements
The results presented thus far were produced by analyzing the full set of observables whenever they were available in terms of illumination, visibility, and eclipse conditions. Nevertheless, it is interesting to study the behavior of the 3σ covariance envelopes when different combinations of orbit data types are considered. For example, Lidar devices are commonly high-power-consuming hardware and are therefore rarely continuously employed for a prolonged time. In the absence of a detailed timeline for the mission's navigation data acquisition and navigation strategy, the different performances that could potentially be achieved by the navigation system processing various sets of observables were analyzed. The effect of the absence of each measurement on the achievable RMS of the position and velocity vector component uncertainties was analyzed. The mean values of the oscillating RMS values and their standard deviations were considered to determine the best achievable performance for each configuration. A time window of a few hours was used to extrapolate these data once the uncertainty 3σ envelopes reached a plateau (as shown in Fig. 10). Figure 13 shows two main phenomena. The Lidar, which has a high frequency of acquisition and low noise, is crucial for this mission. Practically, it is also the device that mainly contributes to lowering the uncertainties regarding Phobos's states. The LOS of the camera appears to provide only minor information for characterizing the moon. It can be observed that the uncertainty levels achieved by processing DSN-only observables are comparable to those recovered without the Lidar data.  Other effects should also be considered. The position and velocity vector uncertainties' RMS is on the same order of magnitude for all combinations of observables when orbiting far from Phobos (QSO-H and QSO-M). Instead, when the spacecraft moves closer without being able to process the Lidar measurements, the difference in performance is evident and appears as a spread between the curves. Poor characterization of the Phobos state vector appears to be capable of significantly degrading knowledge about the MMX state vector as well. An incorrect or insufficient characterization of Phobos introduces perturbations in the sigma point propagation that sum to the process noise stochastic acceleration and worsen the achievable envelope convergence. Lastly, there is a common upward trend in the spacecraft's position and velocity vectors, which is confirmed for all configurations of observables. The overall knowledge of the spacecraft's motion worsens once it moves to QSO-Lb and QSO-Lc, and is more affected by an imperfect estimation of Phobos's gravitational field, owing to the slow-to-decay uncertainties of the moments of inertia.
The suggestion that Lidar is the most important data type for the relative navigation around Phobos can also be observed in Fig. 14. Figure 14 shows how exploiting only the DSN data makes the libration angle non-observable because its uncertainty remains constant    to the a priori values throughout the entire duration of the analysis. Figure 14 also shows that moving closer to the moon does not help improve knowledge regarding the moon's libration motion, which is another reason for the previously mentioned degradation of MMX's state vector knowledge. Upon analysis of the other observables, the presence of the camera's LOS data in the estimation chain generally slightly improves the convergence performance. However, owing to the nature of the observable, the improvement is extremely limited. The exclusion of either of the DSN observables has a similar effect on the convergence rate; however, the presence of the range rate is slightly more important because it is processed more frequently than range measurements.

Comparison between orbit geometry
Another objective of the analysis was to understand whether orbital geometry could influence the characterization of the moon's gravitational field once the mean altitude from Phobos has been selected. The spacecraft's out-of-plane motion of the 3D QSOs is expected to allow a better understanding of the moon's gravitational field [27]. The same should be expected for the swing because of its time-varying altitude profile compared to the nominal QSO solutions. Therefore, we examined how gravitational field knowledge could be refined using different orbit geometries in the lower-altitude orbits. Periodic, 3D, and swing QSO-Lb orbits were used in the second part of the analysis to observe the evolution of the gravitational harmonics coefficients' 3σ envelopes. The resulting uncertainties with respect to the magnitude of the non-normalized harmonic coefficients are shown in Fig. 16.
It is interesting to observe that the periodic QSO is the slowest to converge towards lower values of the gravitational harmonic coefficient uncertainties. This behavior has been observed previously [27] and is confirmed in this analysis, although the difference is not as marked as that previously reported in other studies. A combination of the process noise and the moon's libration motion can partially mask the irregular gravitational field effect on MMX motion. Nevertheless, the out-of-plane motion of the 3D QSO-Lb enables a better estimation of the C 22 gravitational harmonic coefficient, which is better characterized than the other two equatorial orbits by an order of magnitude. The 3D motion allows for a better estimation of Phobos's I zz moment of inertia, which immediately translates into a lower C 22 uncertainty, as shown in Eq. (28). It is believed that the variable altitude characterizing the swing QSO's motion in Phobos's equatorial plane is responsible for the improved performance in observing the I xx and I yy moments of inertia and, therefore, the C 20 coefficient. There is evidently a limit to its observability, but the C 22 uncertainty estimation can still be improved by exposing the spacecraft to different components of Phobos's gravitational acceleration. These considerations should be considered by the flight operation team to maximize the scientific return of the mission. Figure 17 shows that 3D-QSO appears to be the best orbit type for observing the libration angle, albeit by a small amount. The related uncertainty reaches a lower value, and the oscillations of the uncertainty appear flattened with respect to the other orbit geometries. This behavior is again expected because of the out-of-plane motion of the spacecraft, which enables the study of Phobos's libration dynamics from different latitudes. In each result presented so far, it is evident that the observation of every state vector component is fundamental for deeper characterization of the other components. Nevertheless, every element has a different rate of convergence; after the initial fast decay of all the uncertainties' 3σ envelopes, the rate of improvement is dictated by Phobos's MOIs, which are the slowest to decay. Therefore, the decision about which orbit to adopt during the close proximity operation phase should consider the capability of different orbits to characterize the moon's different moments of inertia.

Conclusions
A covariance analysis was performed to investigate the effects of flying quasi-periodic orbits in the vicinity of the Martian moon Phobos. For this analysis, a new dynamic model describing the coupling between the reciprocal gravitational and rotational motion of Mars and Phobos and their effect on MMX's motion was designed. The five periodic quasi-satellite orbits that have currently been baselined for the mission have been transitioned into this mid-fidelity model, along with 3D and swing quasi-periodic orbits. This analysis was performed by comparing the spacecraft's state vector, Phobos's state vector, and the moon's gravitational harmonic coefficient uncertainties' rates of convergence. The developed tool elaborates on the simulated DSN range and range rate measurements, spacecraft onboard Lidar, and LOS observables.
It was shown how the moon's orbital environment could be observed from the five periodic QSOs during the moon's proximity phase. As expected, the more closely the spacecraft orbits the moon, the better the gravitational field can be observed. Nevertheless, there seems to be a lower bound for the achievable estimation of the moon's equatorial moments of inertia when orbiting in the moon's equatorial plane.
Different combinations of orbit data types were also analyzed and considered in the covariance analysis. The spacecraft Lidar was identified as a key component of the mission, being the main source of data for the direct observation of Phobos. As expected, the camera's LOS appeared to be only marginally helpful, whereas the DSNrelated observables (range and range rates) significantly influenced the decay of uncertainties over time.
The orbital geometry was thereafter considered and shown to be important for estimating the gravitational field of Phobos. In particular, the 3D quasi-satellite orbits appear to be more performant in observing the moon's moments of inertia, and should therefore be considered for lower altitudes.
Future developments based on this work may include adding an analysis of navigation techniques based on optical measurements. Phobos's orbital motion around Mars has been characterized quite well over the years. Therefore, autonomous relative navigation strategies could be adopted once the moon's libration motion has been properly characterized and a detailed topographical map is available. Feature-based navigation may help to improve knowledge about the moon's gravitational field, pushing the harmonic coefficients' envelopes plateaus that were observed in this analysis further down.