# Mothership-Cubesat Radioscience for Phobos Geodesy and Autonomous Navigation

Hongru Chen, Nicolas Rambaux, Valéry Lainey, Daniel Hestroffer

► To cite this version:

Hongru Chen, Nicolas Rambaux, Valéry Lainey, Daniel Hestroffer. Mothership-Cubesat Radioscience for Phobos Geodesy and Autonomous Navigation. Remote Sensing, 2022, 14 (7), 10.3390/rs14071619. [https://hal.science/hal-03729621v1](https://hal.science/hal-03729621v1)

HAL Id: hal-03729621

[https://hal.science/hal-03729621v1](https://hal.science/hal-03729621v1)

Submitted on 21 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Article

# Mothership-Cubesat Radioscience for Phobos Geodesy and Autonomous Navigation

Hongru Chen<sup>1,2,3,*</sup> (ID), Nicolas Rambaux<sup>2,3</sup>, Valéry Lainey<sup>2,3</sup> and Daniel Hestroffer<sup>2,3</sup>

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Chen, H.; Rambaux, N.; Lainey, V.; Hestroffer, D. Mothership-Cubesat Radioscience for Phobos Geodesy and Autonomous Navigation. *Remote Sens.* **2022**, *14*, 1619. https://doi.org/10.3390/rs14071619

Academic Editors: Giancarlo Bellucci and Simone Pirrotta

Received: 3 March 2022

Accepted: 21 March 2022

Published: 28 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The knowledge of the interior structure (e.g., homogeneous, porous, or fractured) of Martian moons will lead to a better understanding of their formation as well as the early solar system. One approach to inferring the interior structure is via geodetic characteristics, such as gravity field and libration. Geodetic parameters can be derived from radiometric tracking measurements. A feasible mothership-CubeSat mission is proposed in this study with following purposes, (1) performing inter-sat Doppler measurements, (2) improving the understanding of Phobos as well as the dynamic model, (3) securing the mothership as well as the primary mission, and (4) supporting autonomous navigation, given the long distance between the Earth and Mars. This study analyzes budgets of volume, mass, power, deployment $\Delta v$, and link, and the Doppler measurement noise of the system, and gives a feasible design for the CubeSat. The accuracy of orbit determination and geodesy is revealed via the Monte-Carlo simulation of estimation considering all uncertainties. Under an ephemeris error of the Mars-Phobos system ranging from 0 to 2 km, the autonomous orbit determination delivers an accuracy ranging from 0.2 m to 21 m and 0.05 mm/s to 0.4 cm/s. The geodesy can return 2nd-degree gravity coefficients at an accuracy of 1‰, even in the presence of an ephemeris error of 2 km. The achieved covariance of gravity coefficients and libration amplitude indicates an excellent possibility to distinguish families of interior structures.

**Keywords:** orbit determination; autonomous navigation; planetary geodesy; radioscience; CubeSats; Phobos; small bodies; interior structure

## 1. Introduction

The Martian moons, Phobos and Deimos, are of great scientific interest. The origins of these moons remain unclear. They could possibly be captured C-type asteroids or formed from post-impact debris discs of Mars [1,2]. If they were asteroids, they might play a role in transporting water to terrestrial planets, such as Mars and the Earth [3,4]. Understanding the interior structure of Phobos can put constraints on the formation of Martian moons, as well as that of the early solar system. In addition to the endeavors of Phobos-1 and -2 in the late 1980s, many mission concepts, such as the Phobos-Grunt, PHOTOPRINT, MERLIN, and MMX, for exploring the Martian moons have been proposed and intensively studied by space agencies [4-8]. A mothership-CubeSat radioscience mission is proposed in this paper for probing Phobos and supporting autonomous navigation at the Martian distance from the Earth.

Le Maistre et al. [9] have classified the interior structure of Phobos into four families: rubble pile, heavily fractured, porous, and icy body. The interior structure can be inferred from the moments of inertia (MoI) of Phobos. The MoI is related to the central gravitational parameter, GM, 2nd-degree gravity coefficients, $C_{20}$ and $C_{22}$, and libration amplitude. A number of flybys performed in past missions (e.g., the Viking, Phobos-2, and Mars Express) provided valuable radio tracking data that elucidated the GM of Phobos. One

<footnote n="1">¹ Department of Aeronautics and Astronautics, Kyushu University, 744 Motooka, Fukuoka 819-0395, Japan</footnote>

<footnote n="2">² IMCCE, Observatoire de Paris, CNRS, Université PSL, 77 Av. Denfert-Rochereau, 75014 Paris, France; nicolas.rambaux@obspm.fr (N.R.); lainey@imcce.fr (V.L.); hestro@imcce.fr (D.H.)</footnote>

<footnote n="3">³ IMCCE, Observatoire de Paris, CNRS, Sorbonne Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote>* Correspondence: hongru.chen@aero.kyushu-u.ac.jp</footnote>
latest GM derived radio tracking data during flybys has an uncertainty of 0.3% [10]. The estimation of Phobos' ephemeris and 2nd-degree gravity coefficients done by Jacobson and Lainey [11] shows an 8.7% disagreement on $C_{22}$ from that of the homogeneous model, which corresponds to the shape model derived from observed images [12,13]. While the accuracies of the estimated $C_{20}$ and $C_{22}$ were stated to be 6% and 3% [11], respectively, these results are dependent upon the used quadrupole model of Phobos and the measured libration amplitude (i.e., 1.14° [14]). The libration amplitude estimated based on images has an uncertainty ranging from 2% to 10%, which is in turn dependent on the used ephemeris data [14]. A confident accuracy should be evaluated by considering all uncertainties in the model and measurements.

One latest estimation attempt based on radio tracking data has 1-$\sigma$ uncertainties of 8% and 30% on $C_{20}$ and $C_{22}$, respectively [15]. However, accuracies of 3–5% on $C_{20}$, $C_{22}$, and the libration amplitude are desirable for the identification of the interior structure [9,16]. Nevertheless, the past radio tracking is subject to the flyby geometry relative to the ground station, and the gravity signature is limited by the short duration (i.e., generally <1 h) and far distance (i.e., generally >80 km) of flybys. Consequently, the gravity coefficients $C_{20}$ and $C_{22}$ cannot be resolved at an accuracy sufficient for drawing conclusions on the mass distribution [17]. In addition, Phobos is a secondary body orbiting a primary body, Mars. The gravity and position of the primary can significantly influence the dynamical environment of the secondary as well as the effect of geodesy. Rosenblatt et al. [18] and Yan et al. [19] stated that the error of estimated $C_{20}$ will increase by a few 10% as the ephemeris error tends to a few 100 m. Furthermore, the ephemeris error can go up to a few kilometers by the time of the MMX mission (i.e., 2025–2027) if the ephemeris is not updated. Chen et al. [20] also showed the impact of Phobos' ephemeris errors on spacecraft orbit determination and parameter identification. In summary, ephemeris errors, gravity coefficients, and libration amplitude of Phobos, and orbits of spacecraft are closely correlated in the estimation. It is not only important to enable an observation manner more effective than distant flybys, but also necessary to evaluate the estimation accuracy considering all uncertainties.

Concerning a more effective approach to probing Phobos, there have been several mission concepts and flight plans proposed, the corresponding geodesy performance of which was also evaluated (e.g., [19,21–23]). However, most of the error analyses for the proposed concepts were conducted through analytical covariance analyses. With the high non-linearity related to the flight dynamics and many local optima related to the many variables (e.g., spacecraft orbit states and geodetic parameters) to estimate in the problem, in practice, the theoretical covariance obtained at the global optimum might not be achieved. For the mission proposed in this study, the estimation model is described in Section 4, based on which the theoretical covariance is computed. Massive estimation simulations, which consider all kinds of uncertainties in the model and measurements, are performed to reveal practical covariance (see Section 5). The difference between the theoretical and practical covariances can indicate the functionality of the estimation algorithm and the need for statistical examinations. Despite the difference from the theoretical covariance, the low practical covariance indicates an excellent possibility to distinguish families of interior structures (see Section 6).

The proposed mission concept is designed to feasibly support the desired geodesy and the entire mission. NanoSats have started to play a role in planetary missions. As they are small and low-cost, they can be carried by a mother spacecraft and then deployed to perform audacious flight or landing, or cooperative operations along with the mother spacecraft. There have been several planetary NanoSat missions, such as Philae, MASCOT, MarCo, HERA-Juventas, HERA-Milani [24–26], and ideas (e.g., [23,27–30]). It is proposed in this paper that a CubeSat is carried and deployed by the mother spacecraft aimed for Phobos. The two probes will fly in different quasi-satellite orbits about Phobos, and perform inter-satellite Doppler measurements (see Section 2). The advantage is multiple-fold. First, the measurement is independent of the ground station and thus can support autonomous
navigation, and consequently, enhance the safety of the main mission. Second, the long-duration proximity orbiting will bear an evident signature of the geodetic characteristics of Phobos. As for feasibility, this paper analyzes budgets of volume, mass, power, deployment $Δv$, and link, and the measurement noise of the system, and gives a detailed design of the CubeSat (see Section 3).

## 2. Mission Profile and Orbits

It is proposed that the CubeSat boards the mother spacecraft aimed for Phobos. As precise knowledge on the target is lacking in the beginning, it is preferred that the less indispensable CubeSat is deployed from the mothership to a low-altitude orbit about Phobos, which is sensitive to geodetic parameters of the target, while the mothership stays in a high orbit, which is safe from escape and impact on the target. The two-way Doppler measurement is performed between the two orbiting probes, where the radio signal is triggered by the mothership and reflected by the CubeSat. Section 3.2 will give details on this setup. Based on the inter-satellite measurements, spacecraft orbits can be determined, and the geodetic characteristics of Phobos can be identified. Figure 1 schematically depicts the proposed radioscience mission.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic depiction of the proposed radioscience mission.</figcaption></figure>

The irregularly shaped Phobos has dimensions of 26 × 22 × 18 km. As it is much more smaller than Mars, whose dimension is 6780 km, and orbits closely (i.e., with a semi-major axis of 9377.2 km) around Mars, dynamics in its vicinity is dominated by the gravity of Mars. The sphere of influence of Phobos is small and even below its surface. Therefore, classical Keplerian orbits do not exist around Phobos. The retrograde quasi-satellite orbit (QSO) found in the Mars-Phobos three-body problem is a means for spacecraft to orbit around Phobos in the sense of relative motion [31,32]. Previous work has computed a database of periodic three-dimensional QSO (3D QSO) around Phobos in the circular-restricted three-body problem (CR3BP) [33,34]. Figure 2 presents families of periodic QSO described by the x-amplitude, $A_x$, and z-amplitude, $A_z$. The sensitivity of the orbits to the injection epoch, presumed ground-based orbit determination error (i.e., $1-σ$ 50 m and 3 cm/s on each component), and $Δv$ execution error (i.e., $1-σ$ 1.4 cm/s on each component) has been investigated by propagating the orbits for one week in the presence of the perturbations [34,35]. The color scale in the Figure 2 represents the effective stability of the QSO according to results of sensitivity analyses, which reveals the stability region of bounded orbits around Phobos. The super- and hyper-stable orbits can stay bounded for an acceptable duration (i.e., >7 days) in the high-fidelity model. Low-altitude orbits are generally sensitive to the gravity field of the target, and high-inclination orbits are particularly sensitive to the zonal gravity harmonics [16,22]. Therefore, a low-altitude and high-inclination 29 × 50 × 21 km (in $A_x × A_y × A_z$ format, where $A_y$ is the y-amplitude) 3D QSO, which is marked with a black circle and indicated hyper-stable in Figure 2 and shown to be easily maintainable in Refs. [34,36], is chosen as the candidate science orbit for the CubeSat. For the mothership, a high 100 × 200 km ($A_x × A_y$) planar QSO, which is safe from escape and impact, is adopted.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Solutions of resonant 3D QSO (resonance ratio indicated by *i*: *j*) around Phobos and their effective stability, extracted from Ref. [34].</figcaption></figure>

**3. Budget Analyses and CubeSat Design**

**3.1. Deployment Δv Budget**

Owing to the small capacity of the micro-propulsion system, the Δv cost for deploying the CubeSat from the orbit of the mothership to the low QSO should be low. Ikeda et al. [37] investigated the Δv cost for transfers between QSOs around Phobos and robust injection points in terms of operational safety. They suggested that the injection Δv be executed at leading and trailing points (i.e., at the *y* maximum and minimum in the rotating frame) of the QSOs and the intermediate “swing” QSOs. The swing QSO is the multi-revolution QSO with *A*<sub>*y*</sub> oscillating in the along-track direction of Phobos’ orbit, and connects low and high periodic QSOs. Figure 3 displays an example of a transfer from the 100 × 200 km QSO to the 29 × 46 km QSO in the Phobos-centered rotating frame with Mars on the −*x* axis. Note that the purpose here is simply to have an estimate of the deployment Δv budget. Trajectories and Δv are computed in the simplified CR3BP model. Three impulses are used in this example and all implemented at the leading or trailing points. The transfer is not particularly optimized, while the total Δv is an affordable 16 m/s. To inject the CubeSat to the targeted 3D QSO, an out-of-plane Δv around 3 m/s is required to incline the orbit. Therefore, the deployment requires a Δv approximately at 20 m/s in total. For clear illustration, the shown transfer is the shortest case without revolving around Phobos, which is actually not practical in operation. As Ikeda et al. have noted, since the transfer arcs between QSOs belong to the swing QSOs that keep coming back to the same positions, Δv can be implemented after several revolutions (preferably >2 days) to have a sufficiently long arc for orbit determination at an acceptable operation accuracy (for details, see Refs. [37,38]).
<figure><img src="image_7.png" /><figcaption>Figure 3. Illustration of the deployment from the 100 × 200 km QSO to the 29 × 46 km QSO.</figcaption></figure>

## 3.2. Doppler Measurement and Link Budget

The two-way Doppler setting is selected. In this setting, the signal is generated and emitted by the mothership, then received and transponded by the CubeSat, and finally received and measured by the mothership. The measurement noise consists of the time system noise and the receiver tracking-loop noise. An ultra-stable oscillator (USO) that triggers the radio signals is placed on the mothership. The two-way setting not only simplifies the CubeSat design, but also leads to less time system noise than the one-way setting (i.e., where a less accurate clock is placed on the CubeSat). The thermal noise is related to the received carrier-to-noise density ratio of the radio link, $C/N_0$. To compute $C/N_0$, the link budget should be analyzed. The UHF radio is adopted, which will be explained in the next subsection. The link budget is presented in Table 1 (for the steps computing the link budget, refer to Ref. [39]). Other parameters of the Doppler measurement are presented in Table 2, where the noise in the PLL tracking loop, $\sigma_{f_{PLL}}$ [Hz], is computed from (for details, see Ref. [40]).

$$ \sigma_{f_{PLL}} = \frac{\sqrt{2}}{2\pi t_p} \sqrt{\frac{B_n}{C/N_0} \left[ 1 + \frac{1}{2t_p C/N_0} \right]} \quad (1) $$

As summarized in Table 2, the overall Doppler measurement (or range-rate) noise, $\sigma_{v_r}$, is around 0.09 mm/s.

<table><caption>Table 1. Link budget.</caption><thead><tr><th>Item</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Frequency</td><td><i>f</i></td><td>MHz</td><td>435</td></tr><tr><td>Carrier wave length</td><td>λ<sub>L</sub></td><td>m</td><td>0.6892</td></tr><tr><td>Transmitter output power</td><td><i>P</i><sub>t</sub></td><td>dBW</td><td>-3 <sup>a</sup></td></tr><tr><td>Transmit antenna gain</td><td><i>G</i><sub>t</sub></td><td>dB</td><td>0 <sup>b</sup></td></tr><tr><td>Equiv. isotropic radiated power</td><td><i>EIRP</i></td><td>dBW</td><td>-3</td></tr><tr><td>Propagation path length</td><td><i>D</i></td><td>km</td><td>250 <sup>c</sup></td></tr><tr><td>Space loss</td><td><i>L</i><sub>s</sub></td><td>dB</td><td>-133.18</td></tr><tr><td>Receive antenna gain</td><td><i>G</i><sub>r</sub></td><td>dBi</td><td>0</td></tr><tr><td>Received power</td><td><i>C</i></td><td>dB</td><td>-136.18</td></tr><tr><td>System noise temperature</td><td><i>T</i><sub>s</sub></td><td>K</td><td>635</td></tr><tr><td>Carrier to noise density ratio</td><td><i>C</i>/<i>N</i><sub>0</sub></td><td>MHz</td><td>2.84</td></tr></tbody></table>

<footnote><sup>a</sup> ISIS V/U trceiver. <sup>b</sup> Endurosat UHF antenna. <sup>c</sup> Maximum mothership-CubeSat distance.</footnote>
<table><caption>**Table 2.** Doppler noise budget.</caption><thead><tr><th>Item</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Clock stability</td><td>σ<sub>f<sub>time</sub></sub> / f</td><td>-</td><td>5 × 10<sup>-13</sup> <sup>a</sup></td></tr><tr><td>Code loop noise bandwidth</td><td>B<sub>n</sub></td><td>Hz</td><td>20</td></tr><tr><td>Predetection integration time</td><td>t<sub>p</sub></td><td>s</td><td>5</td></tr><tr><td>PLL thermal noise</td><td>σ<sub>f<sub>PLL</sub></sub> / f</td><td>-</td><td>2.74 × 10<sup>-13</sup></td></tr><tr><td>Overall measurement noise <sup>b</sup></td><td>σ<sub>ν<sub>r</sub></sub></td><td>mm/s</td><td>8.55 × 10<sup>-2</sup></td></tr></tbody></table>

### 3.3. System Sizing

The mission requires constant measurements of Doppler shifts between the mothership and CubeSat via the communication link. However, the limited solar panel area of the CubeSat and the only 51% power generation efficiency, compared to that at 1 AU, constrain the power usage. The UHF radio is adopted in this situation. It allows omni-directional transmission, and thus permits constant Sun pointing for maximum power supply.

Considering these settings, the components and specifications of the CubeSat are decided, as listed in Table 3. All components can be accommodated into a 2U configuration, as shown in Figure 4. The power and mass budgets either based on commercial off-the-shelf specifications or CubeSat development experience are also presented in Table 3. With the largest surface, to which 10 1U solar panels are attached, constantly pointed to the Sun, power generation can cover the consumption during the operating mode. Given a CubeSat dry mass of 2.62 kg, the VACCO 0.25U cold-gas thruster with a total impulse of 93 N·s, specific impulse of 40 s, and wet mass of 676 g can deliver a total Δv of 29 m/s, which covers the need of the deployment, 20 m/s. The total mass of 3.4 kg and volume of 2U should not impose significant loads to the launcher and the mothership. There is still some spare space in the current design. Even 3U and 6U configurations, which are adopted by the CubeSats MarCo [24], Juventas [25], and Milani [26] boarding an interplanetary mothership, are also considered practical. Therefore, depending on the mass and volume budgets allocated to the piggyback, more parts and propulsion capacity may be added to this lowest-cost and most-simplified baseline design, to enhance the robustness and flexibility of the CubeSat mission.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Configuration of the 2U CubeSat.</figcaption></figure>

<footnote><sup>a</sup> General Allen stability of a USO. <sup>b</sup> Root-sum-squared of clock and PLL noises.</footnote>
<table><caption>**Table 3.** Components, mass, and power budgets of the CubeSat.</caption><thead><tr><th>Subsystem</th><th>Components</th><th>Weight, g</th><th>Power, W</th></tr></thead><tbody><tr><td rowspan="2">Structure</td><td>1 2U frame</td><td>390 <sup>a</sup></td><td>-</td></tr><tr><td>Shielding and deployable panels</td><td>523 <sup>a</sup></td><td>-</td></tr><tr><td rowspan="2">Power</td><td>&gt;10 1U solar arrays</td><td>440</td><td>+12.24 <sup>b</sup></td></tr><tr><td>1 battery</td><td>258 <sup>c</sup></td><td>+38.5 h <sup>c</sup></td></tr><tr><td rowspan="5">ADCS</td><td>4 reaction wheels</td><td>220 <sup>a</sup></td><td>-0.8 <sup>a</sup></td></tr><tr><td>1 star tracker</td><td>170 <sup>d</sup></td><td>-1.5 <sup>d</sup></td></tr><tr><td>1 MEMS IMU</td><td>20</td><td>-0.6</td></tr><tr><td>Sun sensors</td><td>- <sup>e</sup></td><td>-0.33</td></tr><tr><td>Communication</td><td>1 UHF antenna set</td><td>85 <sup>f</sup></td><td>-1 <sup>f</sup></td></tr><tr><td></td><td>1 UHF transceiver</td><td>75 <sup>g</sup></td><td>-4 <sup>g</sup></td></tr><tr><td>C&amp;DH</td><td>1 on-board computer</td><td>100</td><td>-0.4</td></tr><tr><td>Propulsion</td><td>1 cold-gas thruster</td><td>676 <sup>h</sup></td><td>-0.25 <sup>h</sup></td></tr><tr><td>Margin</td><td>-</td><td>444 <sup>i</sup></td><td>-1.3 <sup>i</sup></td></tr><tr><td>Total</td><td>-</td><td>3401</td><td>2.0</td></tr></tbody></table>

<sup>a</sup> Based on a lunar CubeSat development experience [41,42]. <sup>b</sup> During the sun-pointing phase with 51% of the efficiency obtained at 1 AU. <sup>c</sup> GOMSpace NanoPower BP4, https://gomspace.com/shop/subsystems/power/ nanopower-bp4.aspx, accessed on 19 February 2019. <sup>d</sup> MAI-SS Space Sextant, https://www.cubesatshop.com/ wp-content/uploads/2016/06/MAI-SS-Specification-10-11-17.pdf, accessed on 1 February 2021. <sup>e</sup> Attached to solar arrays. <sup>f</sup> EnduroSat UHF antenna, https://www.endurosat.com/cubesat-store/all-cubesat-modules/ uhf-antenna/, accessed on 1 February 2021. <sup>g</sup> ISIS V/U tranceiver, https://www.isispace.nl/product/isis-uhf- downlink-vhf-uplink-full-duplex-transceiver/, accessed on 1 February 2021. <sup>h</sup> VACCO End-Mounted Standard MiPS (0.25U) https://www.cubesat-propulsion.com/wp-content/uploads/2015/10/End-mounted-standard- mips.pdf accessed on 1 February 2021. <sup>i</sup> Additional 15% to account for thermal control parts, resistors and coating, and margins of subsystem components.

# 4. Models of Geodesy and Orbit Determination

## 4.1. *Baseline Model*

The dynamic model, as well as the developed orbit propagator, considers ephemerides and gravity of the Sun, Mars, and Phobos. The mar097 ephemeris developed at JPL [43] is adopted as the true ephemeris of Mars and Phobos. The Mars spherical harmonics gravity model is considered up to 10 degrees and orders, whose values are available from the JPL model jgmro_110b. The used Phobos gravity model is up to 3 degrees and orders. The baseline gravity harmonics coefficients of Phobos are computed by the CNES Geodesy Group based on Gaskell's Phobos shape model [44] under the assumption of homogeneous density and GM of $7.1 \times 10^{-4} \text{ km}^3/\text{s}^2$.

The rotational model is computed by Rambaux et al. [45] under the homogeneous assumption, which is expressed as:

$$ \alpha = 317.652 + 1.789 \sin (\mu_1 d + 169.521), \qquad (2) $$

$$ \delta = 52.875 - 1.078 \cos (\mu_1 d + 169.521), \qquad (3) $$

$$ W = 34.781 + \mu_3 d + 1.27 10^{-3} t^2 -1.427 \sin (\mu_1 d + 169.521) - 1.100 \sin (\mu_2 d + 189.271), \qquad (4) $$

where $\alpha$ [°] denotes the right ascension of the figure axis of Phobos, $\delta$ [°] the inclination of its equatorial plane, and $W$ [°] the proper rotation (i.e., $\alpha$-$\delta$-$W$ forms a 3-1-3 rotation sequence with respect to the J2000 ecliptic frame); in addition, $d$ is the integer part of the Julian day, and $t$ is the elapsed seconds on that day. The values of $\mu_1$, $\mu_2$, $\mu_3$ are

$$ \mu_1 = -0.4357344031969911 $$

$$ \mu_2 = 1128.409666972337 $$

$$ \mu_3 = 1128.844884999715 $$
The baseline libration amplitude $\theta$ [°] is -1.1° in Equation (4).

Considering the limited power capacity of the CubeSat, the geodesy campaign should take place during an eclipse-free season. One eclipse-free season is from March 2026 to June 2026 [46]. This period is promising for the Martian mission launched in the year 2024. In the simulation, the initial epoch is set to 2026-4-1 12:00 UTC. As the time of flight of the deployment trajectory is variable, as explained in Section 3.1, and there is no particular requirement on the initial phases of the two orbits, it is assumed that both probes start at the periapsis on the anti-Mars side (i.e., at the x maximum on the x-axis in the rotating frame). Figure 5 presents the evolution of science orbits of the mothership and CubeSat for one week in the described dynamic model. The candidate orbits stay bounded without escape or impact in the high-fidelity model. The Doppler measurement is taken at 1-min intervals. Note that the occultation of Phobos is also considered. In other words, measurements are not valid when Phobos blocks the link. Valid measurements during one week are processed by the estimation routine. In summary, the true orbits of the spacecraft are propagated from true initial orbit states in the high-fidelity dynamic model with the baseline geodetic parameters of Phobos and the mar097 ephemeris, and true range-rates at measurement epochs are generated.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The 100 × 200 km and the 29 × 50 × 21 km spacecraft orbits propagated in the full-dynamic model for 7 days. Orbits are plotted in the Phobos-centered Mars-Phobos rotating frame with Mars on the -x axis.</figcaption></figure>
## 4.2. Error Models

To reveal the influence of ephemeris errors, in the estimation simulation, the mar097 and the NOE-4-2020 developed at IMCCE [47] are adopted as the known ephemerides to express the situations without and with ephemeris errors, respectively. The difference between mar097 and NOE-4-2020 during the period from 2025 to 2027 is around 2 km, mainly along the transverse direction (for more details, see Ref. [47]). The ephemeris error can be approximated by and regarded as a phase shift, $\Delta et$, of 0.9 s.

According to the result obtained by Jacobson and Lainey [11], the non-spherical gravity coefficients of Phobos under the homogeneous assumption can be within a 10% difference from the true situations. The a priori uncertainty of these parameters is set to 10%, and that of GM is set to 3% in this work. Because of the small uncertainty or the small influence, Martian gravity coefficients and parameters other than libration amplitude in Equations (2)–(4) are fixed in the analysis. The a priori uncertainties of the initial orbit states are set to 100 m on each position component and 10 cm/s on each velocity component, which is much greater than the presumed uncertainty of ground-based orbit determination. The Doppler measurement noise $\sigma_{v_r}$, as analyzed in Section 3.2, is conservatively set to 0.1 mm/s. In the Monte Carlo simulation presented in Section 5.3, random errors are generated according to these assumed uncertainties. Table 4 summarizes the uncertainties in the model and measurements. Note that a priori covariances between each of the two uncertain variables are set to zero.

<table><caption>**Table 4.** Summary of error items.</caption><thead><tr><th>Item</th><th>a priori Uncertainty</th></tr></thead><tbody><tr><td>Range-rate measurement</td><td>0.1 mm/s</td></tr><tr><td>Spacecraft position</td><td>(100, 100, 100) m</td></tr><tr><td>Spacecraft velocity</td><td>(10, 10, 10) cm/s</td></tr><tr><td>Phobos' C<sub>00</sub></td><td>3%</td></tr><tr><td>Non-spherical CS coefficients</td><td>10%</td></tr><tr><td>Libration amplitude θ</td><td>0.11°</td></tr><tr><td>Ephemeris error (if present)</td><td>|NOE-4-2020 - mar097|</td></tr></tbody></table>

## 4.3. Theoretical Covariance

The general estimation algorithm optimizes the likelihood of the estimate based on the a priori or predicted covariance and the measurement pre-fit residuals. The *a posteriori* covariance of estimated variables can be expressed as:

$$ \mathbf{P}^* = \mathbf{P}_0 - \mathbf{P}^{e_x e_y} \mathbf{P}^{e_y e_y}^{-1} \mathbf{P}^{e_x e_y} \mathbf{T}, \qquad (5) $$

where $\mathbf{P}^{e_x e_y} (\in \mathbb{R}^{n \times N})$ represents the cross-correlation matrix mapping $n$ uncertain variables $\mathbf{X} (\in \mathbb{R}^{n \times 1})$ (i.e., including initial orbit states of mothership and CubeSat, Phobos' geodetic parameters, and ephemeris errors) to the range-rate measurements $\mathbf{Y} = [v_{r1} \dots v_{rN}]^T$ ($\in \mathbb{R}^{N \times 1}$) obtained at a series of $N$ epochs. $\mathbf{P}_0 (\in \mathbb{R}^{n \times n})$ is the a priori covariance matrix, where the diagonal elements are the squares of values presented in Table 4. $\mathbf{P}^{e_y e_y}$ represents the innovation covariance expressed as:

$$ \mathbf{P}^{e_y e_y} = \mathbf{P}^{yy} + \mathbf{R}^v, \qquad (6) $$

where $\mathbf{P}^{yy} (\in \mathbb{R}^{N \times N})$ represents the output covariance matrix and $\mathbf{R}^v (\in \mathbb{R}^{N \times N})$ the covariance of measurement error. $\mathbf{R}^v$ is a diagonal matrix with $\sigma_{v_r}^2$ in the diagonal entries. $\mathbf{P}^{e_x e_y}$ and $\mathbf{P}^{yy}$ can be obtained using unscented transforms (for details, see Ref. [48]).

The estimation accuracies in two situations, one without ephemeris error and another one under ephemeris errors, are investigated. In this analysis, the mapping from the uncertain ephemeris error to $\mathbf{Y}$ is approximated by the mapping from the ephemeris time shift $\Delta et$. The theoretical estimation accuracy $\mathbf{P}^*$ achieved in one week is computed from
Equation (5). The resulting relative uncertainty is presented in Figure 6. Note that, to show the relative uncertainties of the zero parameters, C₁₀, C₁₁, and S₁₁, 10⁻⁵ is adopted to represent their reference values. As is shown, the accuracies of C₂₀, C₂₂, and libration amplitude, θ, under the ephemeris error are as small as 0.04‰, 0.2‰, and 1‰, respectively. Nevertheless, as the estimation problem involves many uncertain variables and is highly-nonlinear, in practice, the numerical algorithm may not achieve this theoretical accuracy.

<figure><img src="image_10.png" /><figcaption>Figure 6. Theoretical and practical estimation accuracies achieved in one week.</figcaption></figure>

**5. Numerical Estimation Simulations**

This section presents the practical accuracies of Phobos geodesy and spacecraft orbit determination. Random errors are generated based on the a priori uncertainty described in Table 4 and added to the baseline values. Considering 29 uncertain variables (i.e., 6 × 2 orbit states, 3 × 3 CS coefficients, and 1 libration amplitude) and the long computation time to run one estimation routine, 29 × 1000000 samples of the variables are generated as the initial guesses for estimation simulations. Another uncertain factor exists, the ephemeris error. Two situations are simulated, one without ephemeris error in the environment, and the other with ephemeris errors present. The statistical pattern of the covariance is elucidated by Monte-Carlo simulations, as presented in Section 5.3.

Regarding the estimation routine, it is usually difficult to recover both orbit states and geodetic parameters simultaneously. Two scenarios, online orbit determination and offline parameter identification, are planned. The online estimation processes the new measurement to update spacecraft orbit states and the C₀₀ of Phobos using the sequential unscented Kalman filter. In the offline scenario, the batch filter is used to smooth previous estimates and recover the rest of uncertain variables (details of estimation routines will be presented in Ref. [49]).

**5.1. Orbit Determination**

This subsection presents the estimation process for one sample of initial guesses, exhibiting the impact of ephemeris error, and the effect of mitigating this impact. The adopted sampled initial guesses and baseline values (i.e., the assumed true values) are listed in Table 5 for comparison. Note that the a priori and initial guesses of C₁₀, C₁₁, and S₁₁ are set to zero. The process of on-board autonomous navigation (or online orbit determination) under no ephemeris error (i.e., the mar097 ephemeris is employed in the estimation routine) is displayed in Figure 7. The estimation process converges and starts to follow new orbit states in one day. Upon convergence, the average root-sum-squared (RSS) position uncertainty is around 0.2 m, and the RSS velocity uncertainty is 0.05 mm/s.
Figure 8 shows the process of online estimation in the presence of ephemeris errors (i.e., the ephemeris NOE-4-2020 is used). The RSS position uncertainty is significantly increased to 81.6 m, and the RSS velocity uncertainty to 1.63 cm/s. Considering the influence of ephemeris error on the estimation, an ephemeris time error Δet is taken as a consider parameter [49]. Figure 9 shows the estimation process considering the ephemeris time error. The RSS position uncertainty is substantially improved to 20.7 m, and the RSS velocity uncertainty to 0.4 cm/s.

<table><caption>**Table 5.** Presumed truth values and one sample of initial knowledge.</caption><thead><tr><th>Item</th><th>Truth</th><th>Initial Knowledge</th></tr></thead><tbody><tr><td>x<sub>c</sub> [km]</td><td>-16.361</td><td>-16.387</td></tr><tr><td>y<sub>c</sub> [km]</td><td>-10.789</td><td>-10.745</td></tr><tr><td>z<sub>c</sub> [km]</td><td>30.088</td><td>30.128</td></tr><tr><td>v<sub>cx</sub> [m/s]</td><td>1.213</td><td>1.088</td></tr><tr><td>v<sub>cy</sub> [m/s]</td><td>7.249</td><td>7.154</td></tr><tr><td>v<sub>cz</sub> [m/s]</td><td>3.119</td><td>3.045</td></tr><tr><td>x<sub>m</sub> [km]</td><td>-88.829</td><td>-88.880</td></tr><tr><td>y<sub>m</sub> [km]</td><td>-6.438</td><td>-6.470</td></tr><tr><td>z<sub>m</sub> [km]</td><td>46.460</td><td>46.462</td></tr><tr><td>v<sub>mx</sub> [m/s]</td><td>3.592</td><td>3.289</td></tr><tr><td>v<sub>my</sub> [m/s]</td><td>20.561</td><td>20.515</td></tr><tr><td>v<sub>mz</sub> [m/s]</td><td>9.009</td><td>9.133</td></tr><tr><td>C&#773;<sub>00</sub></td><td>1.00000</td><td>0.995321</td></tr><tr><td>C&#773;<sub>20</sub>*</td><td>-0.04757</td><td>-0.05173</td></tr><tr><td>C&#773;<sub>21</sub></td><td>0.00127</td><td>0.00123</td></tr><tr><td>C&#773;<sub>22</sub></td><td>0.02467</td><td>0.02863</td></tr><tr><td>S&#773;<sub>21</sub></td><td>0.00014</td><td>0.00014</td></tr><tr><td>S&#773;<sub>22</sub></td><td>0.00032</td><td>0.00029</td></tr><tr><td>C&#773;<sub>30</sub></td><td>0.00303</td><td>0.00288</td></tr><tr><td>C&#773;<sub>31</sub></td><td>-0.00452</td><td>-0.00464</td></tr><tr><td>C&#773;<sub>32</sub></td><td>-0.00902</td><td>-0.00923</td></tr><tr><td>C&#773;<sub>33</sub></td><td>0.00162</td><td>0.00154</td></tr><tr><td>S&#773;<sub>31</sub></td><td>0.00216</td><td>0.00249</td></tr><tr><td>S&#773;<sub>32</sub></td><td>0.00075</td><td>0.00084</td></tr><tr><td>S&#773;<sub>33</sub></td><td>-0.01360</td><td>-0.01565</td></tr><tr><td>θ [°]</td><td>-1.1000</td><td>-1.1276</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Online determination of orbits of mothership and CubeSat when the ephemeris error is not present.</figcaption></figure>

<footnote>* CS coefficients are normalized values.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Online determination of orbits of mothership and CubeSat when the ephemeris error (i.e., around 2 km) exists and is not considered.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Online determination of orbits of mothership and CubeSat when the ephemeris error (i.e., around 2 km) exists and the ephemeris time error is considered in the estimation.</figcaption></figure>

## Discussion on Orbiting Operations

The obtained accuracy is substantially better than that of the presumed ground-based orbit determination. Moreover, as the convergence time is 1 day, the proposed mission can also support rapid orbit changes independent from ground stations (i.e., without Earth pointing operations nor waiting tens of minutes for commands to travel). As a reference, the permitted orbit maneuver interval in the MMX mission is 2 days, which is associated with the ground-based orbit determination method [38]. Therefore, the proposed orbit determination approach enjoys advances in accuracy, speed, and autonomy.

In addition, for certain mission requirements, the mothership and CubeSat can switch their orbits. For instance, the CubeSat can be first deployed in the low QSO to perform Phobos geodesy, and then deployed to a high QSO serving as a beacon to support the proximity orbiting operation performed by the mothership in a low QSO thereafter. In this case, a bigger cold-gas thruster is required. As reasoned in Section 3.3, slightly bigger configurations are acceptable from a practical perspective. For example, a 0.5U thruster in a 2U CubeSat configuration can deliver a total $\Delta v$ of 49 m/s, and a 0.8U thruster in a 3U CubeSat configuration can deliver a total $\Delta v$ of 68 m/s, which can meet the additional requirement of deployment to a high QSO.
## 5.2. Parameter Identification

The errors of gravity coefficients estimated in the three situations are shown in Figure 10. It can be seen that the ephemeris error significantly interrupts parameter identification. In particular, the libration amplitude $\theta$ cannot be estimated if the ephemeris error is present and not considered in the estimation. With the ephemeris error considered in the estimation, its impact on parameter identification is significantly reduced.

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Errors of gravity coefficients estimated in the three situations; namely, the situation that the ephemeris error is not present, the situation that the ephemeris error exists and is not considered in the estimation, and the situation that the ephemeris error exists and the ephemeris time error is considered in the estimation.</figcaption></figure>

## 5.3. Statistical Covariance

The Monte-Carlo simulation is performed to derive a statistical conclusion. Two situations are considered. One situation is without ephemeris error. The other is under ephemeris errors, and the ephemeris error is considered in the estimation. For each situation, 290 simulations starting with the 290 erroneous samples of uncertain variables are run. Statistical covariance of the gravity coefficients is concluded from differences of the estimates from the true values. The obtained standard deviation is presented in Table 6. The impact of ephemeris errors on geodesy is clearly observed from the differences. Nevertheless, even in the presence of the ephemeris error of 2 km, by considering an ephemeris time error in the estimation, the uncertainty of $C_{00}$ can be kept down to 0.08‰, which is better than the latest implications (e.g., 0.3%, [10]; 0.04%, [15]) from flyby tracking data, suggesting an advanced effectiveness of the proposed dual-sat radiometric observation. The influence of ephemeris error on 2nd-degree gravity identification is around 0.6‰ for $C_{20}$ and 1.7‰ for $C_{22}$. The 1-$\sigma$ uncertainties of $C_{20}$ and $C_{22}$ are hundreds of times smaller than that derived from past flyby tracking data (i.e., 8% and 30%, [15]). The libration amplitude $\theta$, however, becomes almost unobservable under the ephemeris error, as the obtained accuracy is just slightly better than the a priori, 10%. Errors of $C_{20}$ and $C_{22}$ are plotted in Figure 11, along with the concluded covariance ellipses.

<table><caption>**Table 6.** Standard deviation of estimated gravity coefficients elucidated by the Monte-Carlo simulation.</caption><thead><tr><td>Item</td><td>Absence of Ephemeris Errors</td><td>Presence of Ephemeris Errors</td></tr></thead><tbody><tr><td>$C_{00}$</td><td>2.5 × 10<sup>-7</sup></td><td>0.08‰</td></tr><tr><td>$C_{20}$</td><td>0.04‰</td><td>0.65‰</td></tr><tr><td>$C_{22}$</td><td>0.15‰</td><td>1.88‰</td></tr><tr><td>$\theta$</td><td>0.59‰</td><td>8.32%</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Figure 11. Errors and approximated 1-σ and 3-σ certainty ellipses of estimated C20 and C22 under no ephemeris error (blue) and under the ephemeris error (red).</figcaption></figure>

Uncertainties of recovered variables are added to Figure 6 for a clear comparison with the theoretical uncertainties. It can be seen that, the trend of the practical uncertainties generally follows that of the theoretical. While the practical uncertainties are generally greater than the theoretical, the differences are modest. This certifies that estimation routine functions correctly. The difference between the practical and theoretical uncertainties can be related to the non-linearity of the dynamics and the large number of estimated variables. The difference is greater in the situation with ephemeris errors, owing to the unmodeled ephemeris errors in cross-track and radial directions.

## 6. Effect of Inferring Interior Structure

C20 and C22 are related to the MoI of the target. Let Ixx, Iyy, and Izz represent the principle MoI, M the mass, and R0 the reference radius of Phobos. Then, let A, B, and C represent the normalized MoI, Ixx/(MR02), Iyy/(MR02), and Izz/(MR02), respectively. The relationship between 2nd-degree gravity coefficients and MoI is expressed as:

$$C_{22} = (B - A)/4, \qquad (7)$$

$$C_{20} + 2C_{22} = B - C. \qquad (8)$$

The libration amplitude θ [rad] is in related to the flatting γ = (B - A)/C, which is expressed as [50]:

$$\theta = \frac{2e}{1 - \frac{1}{3\gamma}}, \qquad (9)$$

where e = 0.01511 is the orbital eccentricity of Phobos. Equations (7)–(9) can be rewritten as:

$$A = C_{20} + (10 - 24e/\theta)C_{22}, \qquad (10)$$

$$B = C_{20} + (14 - 24e/\theta)C_{22}, \qquad (11)$$

$$C = (12 - 24e/\theta)C_{22}. \qquad (12)$$

Accordingly, the baseline MoIs are A = 0.3537, B = 0.4174, and C = 0.4919.
The small error of GM is negligible in influencing the MoI. To relate errors of C20, C22, and θ to errors of MoI, differential equations are derived as follows,

$$ \begin{bmatrix} \partial A \\ \partial B \\ \partial C \end{bmatrix} = \begin{bmatrix} 1 & 10 - 24e/\theta & 24eC_{22}/\theta^2 \\ 1 & 14 - 24e/\theta & 24eC_{22}/\theta^2 \\ 0 & 12 - 24e/\theta & 24eC_{22}/\theta^2 \end{bmatrix} \begin{bmatrix} \partial C_{20} \\ \partial C_{22} \\ \partial \theta \end{bmatrix}. \quad (13) $$

The Jacobian matrix in Equation (13) indicates that C22 and θ are three to four times more influential than C00 at the same level of relative error, given the baseline values of the parameters. Let M represent the Jacobian matrix. The covariance of the MoI variables on the left-hand side, P_I, is related to the covariance of the parameters on the right-hand side, P_P as:

$$ P_I = MP_P M^T. \quad (14) $$

The previous subsection has deduced P_P from the statistical result. The covariance of MoI is computed from Equation (14) and depicted in Figure 12. When there is no ephemeris error, the accuracy of MoI is as good as 3.3 × 10⁻⁴ (i.e., 0.9‰). When the ephemeris error is present, the 1-σ uncertainty of each normalized MoI is around 0.04 (i.e., 10%), mainly contributed by the θ error. This resolution of MoI cannot directly distinguish families of the interior structure. However, according to Le Maistre et al. [9], the MoI distributions of heavily-fractured and porous-compressed families behave in an opposite way (for details, see Figure 8 of Ref. [9]). Therefore, the near-proportional error distribution exhibited in Figure 12 can distinguish these two families of interior structures.

<figure><img src="image_16.png" /><figcaption>Figure 12. 1-σ uncertainty ellipsoids of inferred moments of inertia projected on the two-dimensional space.</figcaption></figure>
Le Maistre et al. [9] also computed the distribution of different structure families over C20 and C22. The uncertainty ellipse in Figure 11 is placed on the distribution map from Ref. [9], which is shown in Figure 13. Because of the comparably tiny sizes of the uncertainties obtained in this work, only the 3-σ uncertainty ellipse under ephemeris errors is shown. Errors of the estimated C20 and C22 are constrained in this ellipse for the ephemeris error up to 2 km. Given this comparably tiny margin, we can separate the possibility of heavily fractured, porous compressed, and the families in the middle of the map, and distinguish the disrupted-and-reaccreted family from the icy-surface-concentrated and rubble-pile families, regardless of the true values of the parameters.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Map from Ref. [9] on the distribution of structure families over C20 and C22 and the obtained 3-σ uncertainty ellipse (arbitrarily centered on [0.028, −0.046]) under the ephemeris error.</figcaption></figure>

7. Conclusions

A mothership-CubeSat radioscience mission was proposed to probe Phobos and sup-
port autonomous navigation. The two probes are placed in a high (e.g., 100 × 200 km)
quasi-satellite orbit (QSO) and a low three-dimensional (e.g., 29 × 30 × 21 km) QSO about
Phobos, respectively. The operation plan and system design were given to meet require-
ments of power, deployment Δv, and Doppler measurement. The presented analysis of
geodesy and orbit determination takes into account errors in the measurements and model,
including the ephemeris error. Practically achievable estimation uncertainty was revealed
via the Monte-Carlo simulation. The moderate difference between the practical and theo-
retical uncertainties suggests the functionality of the estimation algorithm and the need for
verification via the statistical numerical analysis.

Evident impacts of ephemeris errors on orbit determination and geodesy were observed.
Nevertheless, these impacts were substantially reduced by considering the ephemeris
error in the estimation. The RSS orbit uncertainties under no ephemeris error are 0.2 m
and 0.05 mm/s. The uncertainties increase to 21 m and 0.4 cm/s when an ephemeris
error of approximately 2 km is present. This accuracy is still better than the presumed
accuracy of ground-based orbit determination (i.e., RSS 86 m and 0.5 cm/s). Thus, the
proposed approach can support autonomous orbit determination and rapid (i.e., <1 day)
orbit maneuvers in the proximity of the target, and can tolerate ephemeris errors up to a
few kilometers.
Even under the ephemeris error, the uncertainties of $C_{20}$ and $C_{22}$ recovered from the proposed mission are as small as 0.7‰ and 1.9‰, respectively. However, the *a posteriori* uncertainty of libration amplitude, 8.3%, is not significantly improved from the a priori, 10%. This uncertainty contributes to an uncertainty of 10% in the moment of inertia of Phobos. As the identification of libration is susceptible to the ephemeris error, for high precision of the libration amplitude and the moment of inertia, it is recommended to update the ephemeris of Phobos before the radioscience campaign. Uncertainties achieved in this work are so far the lowest among values obtained in past radioscience campaigns and future mission concepts verified in high-fidelity simulations. The obtained covariance of geodetic parameters can distinguish the porous-compressed interior structure from the heavily-fractured family, and the disrupted-and-reaccreted one from the icy-surface-concentrated and rubble-pile families, regardless of the true values of the parameters.

**Author Contributions:** H.C. and D.H. conceived the ideas; D.H. and H.C. acquired funding; H.C. developed the methodology and software, and conducted the analyses; N.R. and V.L. provided information, study materials, the rotational model and the ephemeris file, and facilitated the access to computing resources; H.C. prepared the original draft; D.H., N.R. and V.L. contributed to the review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present work was initially sponsored by the PSL ESEP Labex for the BIRDY project and is currently sponsored by the Casio Science Promotion Foundation (No. 39-27) and the Education and Research Center for Mathematical and Data Science, Kyushu University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The orbit selection in Section 2 and the true data listed in Table 5 benefit from the results and provided information of the research (e.g., Refs. [33,34,36]) supported by the CNES Research Action (Reference No.: R-S18/BS-0005-039), which are subject to the copyright of CNES.

**Acknowledgments:** The authors appreciate Sebastien Le Maistre (Royal Observatory of Belgium) for allowing the usage of his distribution map in Figure 13. This work benefits from the computer clusters *tycho* and *cerfeuil*, financed and managed by Observatoire de Paris and IMCCE. The first author appreciates the constant support from the informatics department of IMCCE delivered by Stephane Vaillant and Mickaël Gastineau, and also wants to thank Zhuoxi Huo (QXSLab) for an invited visit, Jinsong Ping (NAOC/CAS) for his advice on the Doppler error analysis, and Long Long (CAST) for making the 3D model of the CubeSat.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Abbreviations**

The following abbreviations are used in this manuscript:

3D Three-dimensional

CR3BP Circular-restricted three-body problem

MoI Moment of inertia

QSO Quasi-satellite orbit

RSS Root-sum-squared

**References**

1. Hunten, D.M. Capture of Phobos and Deimos by photoatmospheric drag. *Icarus* **1979**, *37*, 113–123. [CrossRef]

2. Rosenblatt, P.; Charnoz, S.; Dunseath, K.M.; Terao-Dunseath, M.; Trinh, A.; Hyodo, R.; Genda, H.; Toupin, S. Accretion of Phobos and Deimos in an extended debris disc stirred by transient moons. *Nat. Geosci.* **2016**, *9*, 581–583. [CrossRef]

3. Fujimoto, M. MMX (Phobos/Deimos sample return). In Proceedings of the 16th Meeting of the NASA Small Bodies Assessment Group, MD, USA, 12–14 June 2017.

4. Kawakatsu, Y. Mission Design of Martian Moons Exploration (MMX). In Proceedings of the 69th International Astronautical Congress, Bremen, Germany, 1–5 October 2018.

5. Sagdeev, R.Z.; Zakharov, A.V. Brief history of the Phobos mission. *Nature* **1989**, *341*, 581–585. [CrossRef]
6. Marov, M.; Avduevsky, V.; Akim, E.; Eneev, T.; Kremnev, R.; Kulikov, S.; Pichkhadze, K.; Popov, G.; Rogovsky, G. Phobos-Grunt: Russian sample return mission. *Adv. Space Res.* **2004**, *33*, 2276–2280. [CrossRef]

7. Cacciatore, F.; Martín, J. Mission analysis and trajectory GNC for Phobos proximity phase of photoprint mission. *Adv. Astronaut. Sci.* **2015**, *153*, 1321–1340.

8. Murchie, S.; Eng, D.; Chabot, N.; Guo, Y.; Arvidson, R.; Yen, A.; Trebi-Ollennu, A.; Seelos, F.; Adams, E.; Fountain, G.; MERLIN: Mars-Moon Exploration, Reconnaissance and Landed Investigation. *Acta Astronaut.* **2014**, *93*, 475–482. [CrossRef]

9. Le Maistre, S.; Rivoldini, A.; Rosenblatt, P. Signature of Phobos' interior structure in its gravity field and libration. *Icarus* **2019**, *321*, 272–290. [CrossRef]

10. Andert, T.P.; Rosenblatt, P.; Pätzold, M.; Häusler, B.; Dehant, V.; Tyler, G.L.; Marty, J.C. Precise mass determination and the nature of Phobos. *Geophys. Res. Lett.* **2010**, *37*, L09202. [CrossRef]

11. Jacobson, R.; Lainey, V. Martian satellite orbits and ephemerides. *Planet. Space Sci.* **2014**, *102*, 35–44. [CrossRef]

12. Willner, K.; Oberst, J.; Hussmann, H.; Giese, B.; Hoffmann, H.; Matz, K.D.; Roatsch, T.; Duxbury, T. Phobos control point network, rotation, and shape. *Earth Planet. Sci. Lett.* **2010**, *294*, 541–546. [CrossRef]

13. Willner, K.; Shi, X.; Oberst, J. Phobos' shape and topography models. *Planet. Space Sci.* **2014**, *102*, 51–59. [CrossRef]

14. Burmeister, S.; Willner, K.; Schmidt, V.; Oberst, J. Determination of Phobos' rotational parameters by an inertial frame bundle block adjustment. *J. Geod.* **2018**, *92*, 963–973. [CrossRef]

15. Yang, X.; Yan, J.G.; Andert, T.; Ye, M.; Pätzold, M.; Hahn, M.; Jin, W.T.; Li, F.; Barriot, J.P. The second-degree gravity coefficients of Phobos from two Mars Express flybys. *Mon. Not. R. Astron. Soc.* **2019**, *490*, 2007–2012. [CrossRef]

16. Matsumoto, K.; Hirata, N.; Ikeda, H.; Kouyama, T.; Senshu, H.; Yamamoto, K.; Noda, H.; Miyamoto, H.; Araya, A.; Araki, H.; et al. MMX geodesy investigations: Science requirements and observation strategy. *Earth Planets Space* **2021**, *73*, 226. [CrossRef]

17. Pätzold, M.; Andert, T.; Jacobson, R.; Rosenblatt, P.; Dehant, V. Phobos: Observed bulk properties. *Planet. Space Sci.* **2014**, *102*, 86–94. [CrossRef]

18. Rosenblatt, P.; Marty, J.; Dehant, V.; Gurvits, L.; The JIVE Team. *Utilisation de GINS avec les Données de Survol Très Rapproché de Phobos par Mars Express*; In Proceedings of the Journée GINS, Toulouse, France, 12 June 2014.

19. Yan, J.G.; Yang, X.; Ye, M.; Andert, T.; Jin, W.T.; Li, F.; Jin, S.G.; Barriot, J.P. Assessment of Phobos gravity field determination from both near polar and near equatorial orbital flyby data. *Mon. Not. R. Astron. Soc.* **2018**, *481*, 4361–4371. [CrossRef]

20. Chen, H.; Rambaux, N.; Lainey, V.; Hestroffer, D. Mothercraft-cubesat radio measurement for phobos survey. *Adv. Astronaut. Sci.* **2020**, *173*, 387–403.

21. Le Maistre, S.; Rosenblatt, P.; Rambaux, N.; Castillo-Rogez, J.C.; Dehant, V.; Marty, J.-C. Phobos interior from librations determination using Doppler and star tracker measurements. *Planet. Space Sci.* **2013**, *85*, 106–122. [CrossRef]

22. Marty, J.C.; Laurent-Varin, J. Preliminary gravity studies for a Phobos QSO mission. In Proceedings of the 42nd COSPAR Scientific Assembly, Pasadena, CA, USA, 14–22 July 2018; Abstract PSD.1-39-18.

23. Chen, H.; Tomat, L.; Hestroffer, D. CubeSat Complemented Radio Science Small Body Exploration: Survey of Phobos and HO3. In Proceedings of the Interplanetary CubeSat Workshop, Milan, Italy, 28–29 May 2019; p. A.3.4.

24. Klesh, A.T. Inspire and Marco—Technology Development for the First Deep Space Cubesats. In Proceedings of the 41st COSPAR Scientific Assembly, Istanbul, Turkey, 30 July–7 August 2014; Abstract B0.1-18-16.

25. Karatekin, Ö.; Le Bras, E.; Van wal, S.; Herique, A.; Tortora, P.; Ritter, B.; Scoubeau, M.; Moreno, V.M. Juventas Cubesat for the Hera mission. In Proceedings of the European Planetary Science Congress, Online, 13–24 September 2021; p. EPSC2021-750. [CrossRef]

26. Ferrari, F.; Franzese, V.; Pugliatti, M.; Giordano, C.; Topputo, F. Preliminary mission profile of Hera's Milani CubeSat. *Adv. Space Res.* **2021**, *67*, 2010–2029. [CrossRef]

27. Chen, H.; Liu, J.; Long, L.; Xu, Z.; Meng, Y.; Zhang, H. Lunar far side positioning enabled by a CubeSat system deployed in an Earth-Moon halo orbit. *Adv. Space Res.* **2019**, *64*, 28–41. [CrossRef]

28. Cipriano, A.M.; Dei Tos, D.A.; Topputo, F. Orbit Design for LUMIO: The Lunar Meteoroid Impacts Observer. *Front. Astron. Space Sci.* **2018**, *5*, 29. [CrossRef]

29. Asmar, S.W.; Lazio, J.; Atkinson, D.H.; Bell, D.J.; Border, J.S.; Grudinin, I.S.; Mannucci, A.J.; Paik, M.; Preston, R.A. Future of Planetary Atmospheric, Surface, and Interior Science Using Radio and Laser Links. *Radio Sci.* **2019**, *54*, 365–377. [CrossRef]

30. Andrews, D.; Wahlund, J.E.; Kohout, T.; Penttilä, A. Asteroid Prospection Explorer (APEX) Cubesat For the ESA Hera Mission. In Proceedings of the EPSC-DPS Joint Meeting 2019, Geneva, Switzerland, 15–20 September 2019; Volume 2019, p. EPSC-DPS2019-1287.

31. Kogan, Y. Quasi-satellite orbits and their applications. In Proceedings of the Dresden International Astronautical Federation Congress, Dresden, Germany, 6–12 October 1990.

32. Canalias, E.; Lorda, L.; Laurent-Varin, J. Design of realistic trajectories for the exploration of Phobos. In Proceedings of the 2018 Space Flight Mechanics Meeting, Kissimmee, FL, USA, 8–12 January 2018. [CrossRef]

33. Chen, H.; Canalias, E.; Hestroffer, D.; Hou, X. Stability Analysis of Three-dimensional Quasi-satellite Orbits around Phobos. In Proceedings of the 69th International Astronautical Congress, Bremen, Germany, 1–5 October 2018.

34. Chen, H.; Canalias, E.; Hestroffer, D.; Hou, X. Effective Stability of Quasi-Satellite Orbits in the Spatial Problem for Phobos Exploration. *J. Guid. Control Dyn.* **2020**, *43*, 2309–2320. [CrossRef]
35. Canalias, E.; Lorda, L.; Chen, H.; Ikeda, H. Trajectory Design and Operational Challenges for the Exploration of Phobos. *Adv. Astronaut. Sci.* **2020**, *175*, 1493–1507.

36. Chen, H.; Canalias, E.; Hestroffer, D.; Hou, X. Sensitivity Analysis and Stationkeeping of Three-Dimensional Quasi-Satellite Orbits around Phobos. In Proceedings of the 32nd International Syposium on Space Technology and Science, Fukui, Japan, 15–21 June 2019.

37. Ikeda, H.; Mitani, S.; Mimasu, Y.; Ono, G.; Nigo, K.; Kwakatsu, Y. Orbital Operations Strategy in the Vicinity of Phobos. In Proceedings of the International Symposium on Space Flight Dynamics, Matsuyama, Japan, 3–9 June 2017.

38. Oki, Y.; Ikeda, H.; Nishimura, K.; Nakano, M. Operational Safety Analysis on Quasi-Satellite Orbits for Martian Moon eXploration Mission. In Proceedings of the 33rd International Syposium on Space Technology and Science, Oita, Japan, 26 February–4 March 2022.

39. Larson, W.J.; Wertz, J.R. Communications Architecture. In *Space Mission Analysis and Design*, 3rd ed.; Microcosm Press, El Segundo, CA, USA, 2005; Chapter 13.3, pp. 550–570.

40. Kaplan, E.; Hegarty, C.J. *Understanding GPS/GNSS: Principles and Applications*, 3rd ed.; Artech House, Inc.: Norwood, MA, USA, 2017.

41. Long, L.; Jiangkai, L.; Song, H.; Li, L.; Zhu, L.; Chen, H. Preliminary Design and Testing of a Lunar CubeSat System. In Proceedings of the 5th National Symposium on Space Flight Dynamics, Wenchang, China, 27 June–1 July 2017.

42. Yang, C.; Liu, J.; Long, L.; Song, H.; Li, L.; Zhu, L.; Ye, B.; Chen, H.; Zhang, H. Preliminary design and test of attitude control system for lunar CubeSats. *Zhongguo Kongjian Kexue Jishu/Chin. Space Sci. Technol.* **2019**, *39*, 28–35. [CrossRef]

43. Jacobson, R.A. The orbits and masses of the martian satellites and the libration of phobos. *Astron. J.* **2010**, *139*, 668–679. [CrossRef]

44. Gaskell, R. Gaskell Phobos Shape Model V1.0. 2020. Available online: https://arcnav.psi.edu/urn:nasa:pds:gaskell.phobos.shape-model:sdata. (accessed on 2 March 2022). [CrossRef]

45. Rambaux, N.; Castillo-Rogez, J.C.; Le Maistre, S.; Rosenblatt, P. Rotational motion of Phobos. *Astron. Astrophys.* **2012**, *548*, A14. [CrossRef]

46. Nakamura, T.; Ikeda, H.; Kouyama, T.; Nakagawa, H.; Kusano, H.; Senshu, H.; Kameda, S.; Matsumoto, K.; Gonzalez-Franquesa, F.; Ozaki, N.; et al. Science operation plan of Phobos and Deimos from the MMX spacecraft. *Earth Planets Space* **2021**, *73*, 227. [CrossRef]

47. Lainey, V.; Pasewaldt, A.; Robert, V.; Rosenblatt, P.; Jaumann, R.; Oberst, J.; Roatsch, T.; Willner, K.; Ziese, R.; Thuillot, W. Mars moon ephemerides after 14 years of Mars Express data. *Astron. Astrophys.* **2021**, *650*, A64. [CrossRef]

48. Crassidis, J.L.; Junkins, J.L. Sequential State Estimation. In *Optimal Estimation of Dynamic Systems*, 2nd ed.; Chapman & Hall/CRC: Boca Raton, FL, USA; London, UK; New York, NY, USA, 2011; Chapter 3.7, pp. 192–196.

49. Chen, H.; Yoshimura, Y.; Ikeda, H. For orbit determination and geodesy in the non-heliocentric small-body mission. In Proceedings of the International Astronautical Congress, Paris, France, 18–22 September 2022.

50. Duxbury, T.C.; Callahan, J.D. Phobos and Deimos control networks. *Icarus* **1989**, *77*, 275–286. [CrossRef]