# Origin of Phobos and Deimos: Orbital Evolution Shortly after Formation from a Potential Dislocation

Ryan Dahoumane, Kévin Baillié, Valéry Lainey

► To cite this version:

Ryan Dahoumane, Kévin Baillié, Valéry Lainey. Origin of Phobos and Deimos: Orbital Evolution Shortly after Formation from a Potential Dislocation. The Planetary Science Journal, 2025, 6, 10.3847/psj/adcab2. http://hal-05202951

HAL Id: hal-05202951

https://hal.science/hal-05202951v1

Submitted on 7 Aug 2025

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 Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# Origin of Phobos and Deimos: Orbital Evolution Shortly after Formation from a Potential Dislocation

Ryan Dahoumane, Kévin Baillié, and Valéry Lainey

LTE, Observatoire de Paris, Université PSL, Sorbonne Université, Université de Lille, LNE, CNRS, 61 Avenue de l'Observatoire, 75014 Paris, France;
ryan.dahoumane@obspm.fr

Received 2024 August 6; revised 2025 March 17; accepted 2025 April 1; published 2025 May 22

## Abstract

This paper deals with the formation and evolution of Mars' moons, Phobos and Deimos, assuming the dislocation of a larger progenitor as the origin of these moons. The study by R. Hyodo et al. argues that under somewhat simplistic modeling, the postdislocation orbits of Phobos and Deimos inevitably collide within 10,000 yr, leading to their mutual annihilation. These findings are based on N-body simulations, accounting for Mars' $J_2$ and $J_4$ gravitational perturbations and mutual perturbations between the moons. In this paper, we challenge these findings by extending their work. We incorporate important perturbations such as solar perturbations, Mars' axial precession and nutation, and its deformation along three axes. We also extend some of the hypotheses made by R. Hyodo et al. concerning the initial distribution of Phobos and Deimos after the dislocation. Our analysis reveals that including these additional perturbations, as well as the possibility of having more than two fragments after the dislocation, does not alter the ultimate fate of Phobos and Deimos. The moons still converge towards collision within comparable timescales, supporting R. Hyodo et al.'s conclusions that the dislocation hypothesis under the dynamical scenario developed by A. Bagheri et al. has, in the best conditions, about a 10% chance of surviving after the first 100,000 yr following their formation.

*Unified Astronomy Thesaurus concepts: Celestial mechanics (211); Martian satellites (1009); N-body simulations (1083); Natural satellite formation (1425); Planetary system formation (1257)*

## 1. Introduction

The origin of Phobos and Deimos is still a mystery yet to be unraveled. Two main scenarios have emerged from the studies carried out to date: direct capture of Phobos and Deimos and in situ formation. The arguments in favor of a direct separate capture of Phobos and Deimos usually rely on their shape and surface composition similarities with D- or T-type asteroids (P. Rosenblatt 2011). However, considering the small size of Phobos and Deimos, a tidal capture seems very unlikely. Indeed, in such a scenario, both Phobos and Deimos orbits would get high eccentricities and inclinations just after their capture, as a capture on near-circular and near-equatorial orbits is highly unlikely. Assuming that the change of orbits was mainly due to the tidal forces (as it is today), the lowering of their eccentricities and inclinations would require tides that are so strong that they do not seem compatible with what is thought to be the composition of Phobos and Deimos (P. Rosenblatt 2011).

Another way to capture the satellites and lower their inclinations and eccentricities would be through dragging forces. Indeed, D. M. Hunten (1979) has shown that a rotating atmosphere of Mars would be able to capture satellites like Phobos and Deimos and make their orbit quickly evolve to their present positions. However, the main problem with this scenario is that it antedates the formations of Phobos and Deimos to ages when remnants of the rotating disk were still surrounding Mars. This scenario would leave no chance for Phobos to have survived until today, considering its rapid

descent due to tidal forces (V. Lainey et al. 2007). This hypothesis, however, could be a valid explanation for Deimos' origin, assuming the two satellites have different origins.

In situ formation scenarios have also been proposed, for instance, by R. A. Craddock (2011), who defends the idea of a giant impact between Mars and a planetesimal. Indeed, this impact would have produced a debris disk from which Phobos and Deimos could have accreted. This scenario is compatible with a rubble-pile structure for Phobos and Deimos, and helps to explain the various Martian craters that would be due to fallout falls from the accretion disk. A. J. Hesselbrock & D. A. Minton (2017) proposed the hypothesis of reaccretion cycles for the formation of Phobos. Indeed, the proto-Phobos would cross the Roche limit, creating a ring of debris, from which 80% of the mass would fall on Mars and the remaining 20% would reaccrete to form a new generation of Phobos, and so on, diminishing slowly but surely the mass of the Mars innermost satellite. However, in a recent paper, G. Madeira et al. (2023) argued that the reaccretion cycles scenario would imply the existence of a ring close to the Roche limit, which should have been detected in all the previous Martian observations.

Another in situ formation scenario proposed by S. F. Singer (2007) describes Phobos and Deimos as the survivors of the fragmentation of a larger common progenitor satellite. Indeed, the capture of a more massive object is dynamically easier thanks to the greater effects of tides. The capture of this large satellite would have modified Mars' angular momentum, reducing its rotational speed to its current level. At the same time, as the progenitor approached the Roche limit, it would have fragmented, as a result of which the tidal forces would have pulled the more massive fragments toward Mars, leaving only Phobos and Deimos.

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
## 2. The Fragmentation Model

### 2.1. Dislocation and Migration

The scenario we mainly focus on in this paper is the one detailed in A. Bagheri et al. (2021), i.e., the formation of Phobos and Deimos by dislocation and subsequent migration of fragments of a common progenitor satellite. By considering the Mars–Phobos and Mars–Deimos systems as separate systems, they have integrated backwards the semimajor axes, eccentricities, and inclinations of the satellites from their present configuration. These integrations are based on a model of mutual tides between Mars and its satellites, under the assumption that the satellites are always locked in 1/1 spin–orbit resonance. In addition to these tidal forces, they added the physical libration of the satellites, which has the effect of increasing the dissipation due to the tides. The authors estimated the quality factor Q for Phobos based on viscosity estimates and granular friction studies of loose aggregates. The k₂ factor was calculated numerically, using a two-layer model consisting of a solid core and an outer porous layer, each as thick as half the satellite’s radius. For Deimos, these two values were estimated by size-scaling the values of Phobos.

By integrating these equations over several billion years, A. Bagheri et al. (2021) showed that Phobos may have been very close to Mars’ synchronous orbit. The possible high eccentricity of Phobos at that time could even have made the satellite pass above the synchronous radius while having a semimajor axis lower to the synchronous value. This result can be explained by the tidal influence of Mars on Phobos, which has the systematic effect of reducing the satellite’s orbit. If the orbits of Phobos and Deimos are taken far enough back in time, they cross near or beyond the synchronous orbit (more specifically between 5.9 and 6.9 R<sub>Mars</sub>). These crossings would have occurred between 1 Gyr ago for the weakest cases (k<sub>2Pho</sub>/Q<sub>Pho</sub> = 1 × 10<sup>-6</sup> and k<sub>2Dei</sub>/Q<sub>Dei</sub> = 2 × 10<sup>-4</sup>) and 2.7 Gyr for the cases in which Phobos and Deimos are more consolidated (k<sub>2Pho</sub>/Q<sub>Pho</sub> = 6 × 10<sup>-7</sup> and k<sub>2Dei</sub>/Q<sub>Dei</sub> = 1.2 × 10<sup>-4</sup>).

These observations therefore suggest a fragmentation scenario in which a progenitor satellite would generate Phobos and Deimos, which would have then migrated to their current positions thanks to the effects described above.

### 2.2. Caveats

R. Hyodo et al. (2022) published a paper questioning the credibility of A. Bagheri et al.’s (2021) conclusions. Using an *N*-body simulation and the postcollision orbital elements of the satellites calculated by A. Bagheri et al. (2021), R. Hyodo et al. (2022) showed that within 10,000 yr, the dynamical evolution of the system almost systematically results in a collision between Phobos and Deimos, leading to the total destruction of both objects.

They performed 1800 *N*-body simulations of the Mars–Phobos–Deimos system using randomized orbital elements for Phobos and Deimos to estimate the probability of collision. The *N*-body simulations are based on the integration of the equations of motion of the two satellites, which considered mutual gravitational perturbations and the effect of Mars flattening by means of *J*₂ and *J*₄ zonal coefficients. In order to have initial conditions unlikely to generate collisions, the authors have considered only cases in which, after fragmentation, the orbits of both satellites have only one point of contact, which corresponds to the apoapsis of Phobos and the periapsis

of Deimos (i.e., $a_{Pho}(1 + e_{Pho}) = a_{Dei}(1 - e_{Dei})$ and $|\omega_{Pho} - \omega_{Dei}| = \pi$). In spite of these considerations, a collision between both fragments happens in more than 94% of the simulations. From impact simulations, the authors demonstrate that such collisions (100–300 m·s⁻¹) would lead to annihilation of both satellites, creating a ring of debris. Based on these conclusions, it seems that the model proposed by A. Bagheri et al. (2021) is unlikely.

### 2.3. Limitations and Further Constraints

#### 2.3.1. Postfragmentation Orbital Elements

Figure 1 from A. Bagheri et al. (2021) as well as the conditions on the apsides described above allowed R. Hyodo et al. (2022) to deduce the following constraints on the postfragmentation orbital elements of Phobos and Deimos:

$$ 
\begin{align*}
a_{Pho} & \in [5R_{Mars}, 6R_{Mars}] & a_{Dei} & \in [6.5R_{Mars}, 7.5R_{Mars}] \\
e_{Pho} & \in [0.15, 0.35] & e_{Dei} & \in [0, 0.2] \\
i_{Pho} & = 0.021 \text{ rad} & i_{Dei} & = 0.015 \text{ rad} \\
\Omega_{Dei} & = \Omega_{Pho} & |\omega_{Dei} - \omega_{Pho}| & = \pi \text{ rad} \\
a_{Pho}(1 + e_{Pho}) & = a_{Dei}(1 - e_{Dei}),
\end{align*}
$$

where *a*, *e*, *i*, *ω*, and *Ω* are, respectively, the semimajor axes, eccentricities, inclinations, arguments of the periapsis, and longitudes of the ascending node. This orbital elements distribution can be seen in Figure 1.

We can make several comments on the choice of these constraints:

1. Figure 1 by A. Bagheri et al. (2021), showing the evolution of the eccentricity of Phobos and Deimos, actually indicates a minimum postfragmentation eccentricity of Phobos of slightly less than 0.15. For their simulations, R. Hyodo et al. (2022) should therefore probably have distributed the eccentricities of Phobos between 0.14 and 0.35 in order to be conservative.

2. It is explicitly stated in A. Bagheri et al. (2021) that the orbits intersect at a distance of 5.9–6.9 Mars radii. This constraint does not seem to have been taken into account by R. Hyodo et al. (2022) in their distribution.

This second point significantly changes the distribution of the initial orbital elements of the two moons, in particular by greatly reducing the mean eccentricity of Deimos. Note that circular orbits increase the chances of collisions, since zero eccentricity virtually eliminates the relative precession of the apsides, and only the precession of the nodes can separate the orbits. From these new orbital distributions presented in Figure 1, we can expect a drop in collision probability, particularly for high values of Phobos’ semimajor axis, when this new constraint drastically reduces the number of circular orbits.

Furthermore, in order to minimize the chances of collision, R. Hyodo et al. (2022) had assumed that the apoapsis of Phobos is equal to the periapsis of Deimos. This is a good idea; however, we will consider a situation that minimizes further the chances of collisions by slightly increasing the distance between satellites in the initial conditions. Instead of making the center of the satellites coincide using the condition $a_{Pho}(1 + e_{Pho}) = a_{Dei}(1 - e_{Dei})$, we will consider the case where the two moons are barely in contact after fragmentation: $a_{Pho}(1 + e_{Pho}) + R_{Pho} = a_{Dei}(1 - e_{Dei}) - R_{Dei}$, with R being
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distribution of the initial semimajor axes and eccentricities of the simulations. Each colored box represents the orbital elements distribution allowed by the constraints mentioned in Section 2.3.1. The “H22” refers to the one used by R. Hyodo et al. (2022) in their simulations, while the “This Paper” shows those used in the present paper. Note that the green color is the superposition of yellow and light blue, while the purple is the overlap of red and dark blue.</figcaption></figure>

the radius of each moon. In this way, we expect the number of collisions at the start of the simulation to be greatly reduced, as any small perturbation would prevent them. This new orbital elements distribution is shown in Figure 1 and is based on the following constraints:

$$a_{\text{Pho}} \in [5R_{\text{Mars}}, 6R_{\text{Mars}}] \quad a_{\text{Dei}} \in [6.5R_{\text{Mars}}, 7.5R_{\text{Mars}}]$$

$$e_{\text{Pho}} \in [0.14, 0.35] \quad e_{\text{Dei}} \in [0, 0.2]$$

$$i_{\text{Pho}} = 0.021 \text{ rad} \qquad i_{\text{Dei}} = 0.015 \text{ rad}$$

$$\Omega_{\text{Dei}} = \Omega_{\text{Pho}} \qquad |\omega_{\text{Dei}} - \omega_{\text{Pho}}| = \pi \text{ rad}$$

$$5.9R_{\text{Mars}} \le a_{\text{Pho}}(1 + e_{\text{Pho}}) + R_{\text{Pho}} = a_{\text{Dei}}(1 - e_{\text{Dei}}) - R_{\text{Dei}} \le 6.9R_{\text{Mars}}.$$

### 2.3.2. Limited Physical Effects

In addition, we propose to increase the complexity of the physical model used by R. Hyodo et al. (2022) in their paper to test the robustness of their results. The authors justified the simplicity of their model by explaining that the most important perturbation apart from the flattening of Mars is ejection, which was said to have no impact on the result. We propose to verify this claim. We extended the R. Hyodo et al. (2022) model by adding the gravitational perturbations due to the Sun, the dynamical triaxiality of Mars (coefficients $c_{22}$ and $s_{22}$), and the precession and nutation of Mars.

Solar perturbations are expected to influence the nodal precession of satellite orbits, particularly when they are orbiting far from Mars. This effect is expected to intensify when increasing Martian obliquity. Additionally, the triaxial shape of Mars may induce resonant interactions with the satellites when they approach spin-orbit resonance. Furthermore, the precession of Mars' rotational axis could induce periodic variations in the satellites' orbital inclinations. However, this influence is likely to be minor, as the precession of the satellites' ascending nodes occurs at a faster rate than the planetary precession (see P. Goldreich 1965 and P. Gurfil et al. 2007).

For the sake of completeness, we also added the $J_3$ flattening term for Mars. This term was neglected by R. Hyodo et al. (2022) because the low inclinations of Phobos and Deimos make the $J_3$ perturbations negligible.

Additionally, we included tidal effects in our modeling. Within such short timescales (tens of thousands of years), the orbital change due to tides within Mars, as well as within

Phobos and Deimos, is expected to be negligible. However, in order to be as exhaustive as possible, we have conducted additional simulations where we considered extreme values for the $k_2/Q$ of Mars, Phobos, and Deimos significantly higher than present values. We neglect the relativistic effect as the relativistic shift of Phobos periapsis is about $3 \cdot 10^{-7} \text{rad} \cdot \text{yr}^{-1}$ based on P. Arrighi & G. Dowek (2015) ($\Delta\omega = \frac{24\pi^2 a^2}{T^2 c^2 (1 - e^2)}$ with $\Delta\omega$, $c$, and $T$, respectively, the shift in periapsis, the speed of light, and the Newtonian orbital period). As a comparison, based on Equation (14) from R. Hyodo et al. (2022), the precession of Phobos periapsis due to Mars oblateness is about $0.7 \text{ rad} \cdot \text{yr}^{-1}$, 5 to 6 orders of magnitude greater.

## 3. Methodology

### 3.1. Dynamical Modeling

#### 3.1.1. Equations of Motion

In a planetocentric frame, the acceleration undergone by each moon considered as point masses can be written as (V. Lainey et al. 2004)

$$\ddot{r}_i = -\frac{G(m_0 + m_i)r_i}{r_i^3} + \sum_{j=1, j \ne i}^{N} Gm_j \left( \frac{\mathbf{r}_j - \mathbf{r}_i}{r_{ij}^3} - \frac{\mathbf{r}_j}{r_j^3} \right) + G(m_0 + m_i) \nabla_i V_{i\hat{0}} \quad (1)$$

The subscripts $i$ and $j$ represent any of the $N$ bodies in the problem, with the exception of the central body, which is symbolized by the subscript 0. $G$ is the universal gravitational constant, $\mathbf{r}_i$ represents the position vector of body $i$ in the reference frame described above, $r_i$ represents the norm of this vector, and $r_{ij}$ represents the distance between body $i$ and body $j$. $\nabla$ designates the gradient operator, and $V_{i\hat{0}}$ represents the gravitational term due to the shape of the central body exerted on the $i$ body. As in V. Lainey et al. (2004), we divided this last term into two parts: $V_{i\hat{0}} = V_{i\hat{0}}^{(1)} + V_{i\hat{0}}^{(2)}$, where $V_{i\hat{0}}^{(1)}$ is the term due to the flattening of the central body and $V_{i\hat{0}}^{(2)}$ is the term due to the sectorial terms.

The derivative of the term due to the flattening $J_2$, $J_3$, and $J_4$ can be written as follows (Equation (1.4) from
V. Lainey et al. 2004):

$$ \nabla_i V_{i\hat{0}}^{(1)} = -\sum_{n=2}^{\infty} \frac{(R)^n J_n}{r_i^{n+2}} \left[ \left( \mathbf{k} - \frac{\mathbf{r}_i}{r_i} \sin \phi_i \right) P_n'(\sin \phi_i) - (n+1) \frac{\mathbf{r}_i}{r_i} P_n(\sin \phi_i) \right]. \quad (2) $$

Expanding the expression to order 4, we find

$$ \nabla_i V_{i\hat{0}}^{(1)} = \frac{(R)^2 J_2}{r_i^4} \left[ \frac{\mathbf{r}_i}{r_i} P_3'(\sin \phi_i) - k P_2'(\sin \phi_i) \right] \\ + \frac{(R)^3 J_3}{r_i^5} \left[ \frac{\mathbf{r}_i}{r_i} P_4'(\sin \phi_i) - k P_3'(\sin \phi_i) \right] \\ + \frac{(R)^4 J_4}{r_i^6} \left[ \frac{\mathbf{r}_i}{r_i} P_5'(\sin \phi_i) - k P_4'(\sin \phi_i) \right], \quad (3) $$

with $k$ the rotation vector of the central body, $R$ the equatorial radius of the central body, $\phi_i$ the latitude of the body $i$ with respect to the equator of the central body, $P_n$ the Legendre polynomial of degree $n$, and $P'_n$ the derivative of the Legendre polynomial of order $n$ with respect to $\sin \phi$. The sine of the latitude can easily be found as follows:

$$ \sin(\phi_i) = \frac{\mathbf{r}_i \cdot \mathbf{k}}{r_i}. \quad (4) $$

As mentioned previously, we decided to include the term related to the $J_3$ coefficient in the calculation of the potential for the sake of completeness. However, even without this term, we have a slightly different expression from that given by R. Hyodo et al. (2022) for three reasons.

1. First, the expression for the potential due to the flattening of the central body is different (Equation (3)). This is because the authors assumed a fixed rotation axis and thus defined a reference frame xyz with xy the Martian equatorial plane and z the rotation axis of the central body. Given that we are going to consider the precession and nutation of the central body's rotation axis, we are not in the same type of frame and therefore have a different, more general expression. Indeed, by developing Equation (3), considering $k = (0, 0, 1)^T$ and neglecting the term due to $J_3$, we end up with the same expression of the potential due to the flattening of the central body as that given by R. Hyodo et al. (2022).

2. Second, they have neglected the mass of the body $i$ when summing it up with the mass of Mars in the first term of the Equation (1). This could be done, as Phobos and

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Representation of the fixed and moving frames. The first frame, which is fixed, is defined by arbitrary axes ($X$, $Y$, $Z$). After three rotations with the angles $(I, \psi, \chi)$, this fixed frame can be turned into the rotating frame ($X'$, $Y'$, $k'$), where $X'$ is toward the first meridian and $k'$ is the perpendicular to the equator. We also represented $\bar{X}$, which is directed toward the ascending equatorial node.</figcaption></figure>

3. Finally, it also appears that the authors omitted the indirect terms in their calculation of the force. Fortunately, this omission does not seem to affect the results since, as can be seen from Figures 3 and 4, the results of their simulations can be found from our formulation of the force, which includes the indirect terms.

Let us now express the derivative of the term depending on $c_{22}$ and $s_{22}$, due to the disparities in longitude of the Mars mass distribution (sectorial terms). This can be written as (V. Lainey et al. 2004):

$$ \nabla_i V_{i\hat{0}}^{(2)} = \frac{(R)^2}{r_i^4} \left[ \frac{\mathbf{r}_i}{r_i} (15 \sin^2 \phi_i - 9) - 6k \sin \phi_i \right] \\ \times [c_{22} \cos 2\lambda_i + s_{22} \sin 2\lambda_i] \\ + \frac{(R)^2}{r_i^3} 6(1 - \sin^2 \phi_i) \frac{x_i' \nabla y_i' - y_i' \nabla x_i'}{x_i' 2 + y_i' 2} \\ \times [-c_{22} \sin 2\lambda_i + s_{22} \cos 2\lambda_i], \quad (5) $$

where $\lambda_i$ is the longitude of the position of the body $i$ and $(x'_i, y'_i, z'_i)$ the coordinates of the body $i$ in the rotating frame centered on the central body ($P_0$, $X'$, $Y'$, $k'$), with $x'_i$ the equatorial coordinate pointing toward the first meridian, $z'_i$ the coordinate pointing along the axis of rotation of the central body, and $y'_i$ completing the direct trihedron (see Figure 2). The transition matrix $M$ between the coordinates of this rotating frame and the fixed frame ($P_0$, $X$, $Y$, $Z$) such that $(x'_i, y'_i, z'_i)^T = M(x_i, y_i, z_i)^T$ is the following (V. Lainey 2002; V. Lainey et al. 2004):

$$ M = \begin{pmatrix} \cos \chi \cos \psi - \sin \chi \sin \psi \cos I & \cos \chi \sin \psi + \sin \chi \cos I \cos \psi & \sin \chi \sin I \\ -\sin \chi \cos \psi - \cos \chi \sin \psi \cos I & -\sin \chi \sin \psi + \cos \chi \cos \psi \cos I & \cos \chi \sin I \\ \sin I \sin \psi & -\cos \psi \sin I & \cos I \end{pmatrix}, \quad (6) $$

Deimos' masses are negligible when compared to Mars. However, as we include the Sun in our modeling, $m_i$ can not be neglected in this paper.

with $I$, $\psi$, and $\chi$ being, respectively, the inclination of the rotation axis (obliquity), the longitude of the ascending node, and the angle between the ascending node and the prime meridian.
### 3.1.2. Computing the Precession and Nutation of Mars

The precession and nutation motions of a body orbiting the Sun are given by P. Gurfil et al. (2007). It is possible to describe the evolution of Mars' axis of rotation as follows:

$$ \frac{dk}{dt} = \alpha(\mathbf{n} \cdot \mathbf{k})(\mathbf{k} \times \mathbf{n}), \qquad (7) $$

where $\mathbf{k}$ denotes the Martian spin axis, $\mathbf{n}$ the normal to the orbital plane, and $\alpha$ a parameter proportional to the oblateness factor $J_2$. The $\mathbf{n}$ vector can be written as

$$ \mathbf{n} = (q, -p, \sqrt{1 - p^2 - q^2})^T, \qquad (8) $$

with

$$ q = \sin I_{\text{orb}} \sin \Omega_{\text{orb}} \quad p = \sin I_{\text{orb}} \cos \Omega_{\text{orb}}, \qquad (9) $$

where $I_{\text{orb}}$ is the inclination of the Martian orbit with respect to a fixed plane and $\Omega_{\text{orb}}$ is the ascending node with respect to this same plane. It is trivial to calculate these values directly from the relative positions and velocities of Mars and the Sun; however, for the sake of precision, similarly to P. Gurfil et al. (2007), we use W. R. Ward's (1979) model, which describes the evolution of these angles as a function of time in the fixed reference frame B1950. We therefore use Equations (20) and (21) from P. Gurfil et al. (2007) to calculate the coefficients $q$ and $p$.

### 3.1.3. Tides in Mars

To take into consideration the effect of Martian tides on Phobos and Deimos, we have used the expression given in Equation (3) from V. Lainey et al. (2007), which gives the force acting on a moon induced by the bulge created by this moon on the planet. Considering that Mars is the center of our referential, we can write

$$ \ddot{\mathbf{r}}_i^T_0 = -\frac{3k_2Gm_i(R)^5}{r_i^8} \times \left( \mathbf{r}_i + \Delta t \left[ \frac{2\mathbf{r}_i(\mathbf{r}_i \cdot \mathbf{v}_i)}{r_i^2} + (\mathbf{r}_i \times \mathbf{\rho} + \mathbf{v}_i) \right] \right), \qquad (10) $$

with $\mathbf{v}_i$, $k_2$ and $\Delta t$, respectively, the velocity of the body $i$, the planet's Love number, and the time shift between the excitation of the planet and its viscoelastic response. $\mathbf{\rho}$ represents the rotation vector of the planet such that

$$ \mathbf{\rho} = k \dot{\mathbf{\chi}}. \qquad (11) $$

Moreover, Equations (1) and (2) from V. Lainey et al. (2007) give us

$$ \Delta t = \frac{T_0^{\text{rot}} T_i \arcsin(1/Q)}{4\pi |T_0^{\text{rot}} - T_i|}, \qquad (12) $$

with $T_0^{\text{rot}}$, $T_i$, and $Q$ the planet's rotation period, the satellite's revolution period around the planet, and the planet's tidal quality factor.

In this work, we will assume $k_2$ and $Q$ as constants, as the point is to check if extreme values can alter the results rather than having an accurate representation of the phenomenon. Moreover, because of the small size of the satellites compared to the planet, we neglect the influence of Phobos and Deimos on the Mars rotation speed $\dot{\mathbf{\chi}}$.

### 3.1.4. Tides in Phobos and Deimos

Finally, the last effect we added in the modeling of our problem is the force exerted by the moons on the planet due to the tidal bulges of the moons. In an analog way to Equation (10), the acceleration of the body $i$ due to this force in a planetocentric referential can be written as

$$ \ddot{\mathbf{r}}_i^T = -\frac{3k_2i Gm_0 (R_i)^5}{r_i^8} \times \left( \mathbf{r}_i + \Delta t_i \left[ \frac{2\mathbf{r}_i(\mathbf{r}_i \cdot \mathbf{v}_i)}{r_i^2} + (\mathbf{r}_i \times \mathbf{\rho}_i + \mathbf{v}_i) \right] \right), \qquad (13) $$

and

$$ \Delta t_i = \frac{T_i^{\text{rot}} T_i \arcsin(1/Q_i)}{4\pi |T_i^{\text{rot}} - T_i|}, \qquad (14) $$

with $k_{2i}$, $Q_i$, $R_i$, $\rho_i$, $T_i^{\text{rot}}$, and $\Delta t_i$, respectively, the Love number, tidal quality factor, equatorial radius, rotation vector, and time lag of the body $i$ (satellite). We still consider $k_{2i}$ and $Q_i$ as constant for the same reasons we did with $k_2$ and $Q$. As we want to maximize the effect of the tides in the satellite, we will consider that, at any given time, the rotation of the satellites is in the same plane as their revolution around the planet.

Unlike the planet, due to their small size, the variation of the rotation speed of the satellites can not be neglected when computing the acceleration due to their tidal deformation. In other words, we cannot consider the rotation speed $\rho_i$ of the satellite as constant. For this reason, we will be simultaneously integrating the position rotation of the satellites. To do so, we will be using the Equation (7) from E. Y. Aleshkina (2009),

$$ -\frac{d\rho_i}{dt} = K_{i,i} \left(\frac{1}{r_i}\right)^6 + K_{i,i} \left(\frac{a_i}{r_i}\right)^3 \sin 2(\theta_i - f_i), \qquad (15) $$

with $a_i$ the semimajor axis of the satellite, $f_i$ its true anomaly, and $\theta_i$ the angle between the axis of the moment of inertia $A$ of the satellite and its eccentricity vector ($\rho_i = \frac{d\theta_i}{dt}$). The first term of the expression is due to the tidal bulge, while the second term is related to the inertia of the satellite. If we considered satellites as point masses, the second term would not appear. However, as Phobos and Deimos are not spherical, it must be taken into account. The expressions for the tidal coefficient $K_{i,i}$ and the inertial coefficient $K_{i,i}$ given in E. Y. Aleshkina (2009) are

$$ K_{i,i} = \frac{15}{4} \frac{Gm_0^2 R_i^3}{m_i Q_i} k_{2i} \qquad (16) $$

$$ K_{i,i} = \frac{3}{2} \frac{B - A}{C} n_i^2, \qquad (17) $$

with $A$, $B$, and $C$ the main moments of inertia of the satellite, and $n_i$ its mean motion. As in E. Y. Aleshkina (2009), we will assume that $C/m_i R_i^2 \approx 0.4$. Moreover, the $c_{22}$ term of the satellite can be written (C. F. Yoder 1995) $c_{22} = \frac{1}{4} \frac{B-A}{m_i R_i^2}$. From
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Cumulative fraction of collisions as a function of simulation time. We made different simulations considering different configurations. For each of these configurations, we made 600 simulations for each value of $a_{Phobos}$. The “Original” curves are a pure reproduction of R. Hyodo et al. (2022) results, and “New OE” curves have the same physical effects as those considered by R. Hyodo et al., but with the new initial orbital elements repartition detailed in this paper. Then, in the “J3” curves, we also added the oblateness term $J_3$, after which we also added the Sun gravitational perturbations in the “J3 + Sun” curves. As mentioned in its name, we then added the precession in the “J3+Sun+Precession” curves, and finally, we added the Martian sectoral terms $C_{22}$ and $S_{22}$ in the “Full” simulation.</figcaption></figure>

this, the inertial coefficient can be rewritten as

$$K_{l,i} \approx 15c_{22}n_i^2. \qquad (18)$$

## 3.2. Numerical Integration: Hybrid Radau Integrator

We now aim to solve numerically the equations of motion described above. In order to combine the speed of Everhart’s original RA15 code (E. Everhart 1985) with the improved accuracy provided by Rein and Spiegel’s implicit integrator with adaptive time stepping, 15th order (IAS15) (H. Rein & D. S. Spiegel 2014), we developed a hybrid integrator that can switch methods between RA15 and IAS15 at each time step depending on the distance between the different bodies within an *N*-body simulation. This means that when the bodies are far apart, the classical Radau is used for fast integration. However, when two bodies are sufficiently close, the integrator used is the IAS15, which guarantees a higher accuracy for close encounters and ensures that no collisions are missed.

To do so, after each time step, a check is made on the distance between each preselected body. If the distance between these objects is less than a threshold $\epsilon$, the integration of the next time step is made by the IAS15, otherwise, the Radau is used. We have set the threshold to activate the IAS15 when the lighter body is susceptible to enter the Hill sphere of the heavier.

Let the *i* subscript symbolize the heavier body and the *j* subscript symbolize the lighter one. Let $v_{j \to i}$ be the velocity vector of *j* relative to *i*, and $r_{hill}(i)$ the radius of the Hill sphere of the body *i*. Assuming both bodies are spherical, let $R_i$ and $R_j$ be the radii of each body and *d* the distance between them. Then, considering that the velocity magnitude is constant during the timestep, we can set

$$\epsilon = r_{hill}(i) + (R_i + R_j) + \|v_{j \to i}\| * (t_{n+1} - t_n). \qquad (19)$$

## 3.3. Detecting Collisions in the Simulation

Collision detection plays a major role in the simulations. The most straightforward way to detect a collision between two

bodies in an *N*-body problem consists in checking the distance *d* between them at each timestep:

$$d \le R_i + R_j \Rightarrow \text{collision.} \qquad (20)$$

This basic method presents, however, a flaw: it would miss a collision that occurs between two successive time steps. Nevertheless, our testing indicates that the IAS15 part of the integrator we are using is very efficient in reducing the timestep during close approach, so we will not need any additional collision detection condition. In addition, we want to avoid false positives as much as possible, as our simulations aim to set a higher bound on the survival rate.

# 4. Results

## 4.1. Changing the Initial Distribution and Extending the Simulation Time

We first want to see how our new orbital elements distribution (described in Section 2.3.1) affects the collision probability while using the same physical process as R. Hyodo et al. (2022). By comparing the blue and the green curves from Figure 3, we notice that, despite notable differences at the beginning of the simulation, we obtain a result very similar to that of R. Hyodo et al. (2022) after 100,000 yr of evolution: more than 94% of the simulations have undergone a collision. However, there is a very sharp drop in the number of collisions in the short term, which can be explained by the fact that we have shifted the positions of the orbits of Phobos and Deimos so that they are separated by the sum of their physical radii. This new separation prevents the two bodies from colliding in their first revolution, but the mutual perturbations of Phobos and Deimos have the secular effect of bringing their orbits closer together, leading to a similar result as the one obtained by R. Hyodo et al. (2022).

As for the effect induced by the new criterion (fragmentation distance between 5.9 and 6.9 Mars radii), which was omitted by R. Hyodo et al, it has a very weak effect on the evolution of the number of collisions, but has a somewhat stronger effect on
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Cumulative fraction of collisions as a function of impact velocity. We did the same plot as in R. Hyodo et al. (2022) to compare their results with our new simulations. These results come from the same simulations as those presented in Figure 3. The plain, dashed, and dotted lines represent represent, respectively, the reproduction of R. Hyodo et al.'s (2022) results, the results from a physical model similar to R. Hyodo et al. (2022) but with new orbital elements (see Figure 1), and the simulations with our new physical modeling and orbital elements distribution. The colors light, middle, and dark gray represent, respectively, the simulations in which we considered the initial semimajor axis of Phobos as equal to 6, 5.5, and 5 Mars radii.</figcaption></figure>

the cases where $a_{Pho} = 6R_{Mars}$. In this case, the R. Hyodo et al. (2022) conditions tend to generate lower Deimos eccentricities than those of the new configuration. Calculating the median Deimos eccentricity for the 6 Mars radii case typically gives a value of around 0.02 for the R. Hyodo et al. (2022) distribution, compared with a value of around 0.04 for the new distribution. It is possible that this reason is partly responsible for the reduction in the amplitude of the first collision peak in the new distribution of orbital elements in the case where $a_{Pho} = 6R_{Mars}$.

Despite this, it seems that a distribution even more unfavorable to collisions than the previous one is not enough to significantly reduce the probability of Phobos and Deimos colliding, as the blue and green curves from Figure 3 are very close to each other after $10^5$ yr.

## 4.2. Adding Physical Effects

We now investigate the impact of each of the physical effects added to the simulation (results available in Figure 3).

As expected, given the very low inclinations of Phobos and Deimos, the effect of $J_3$ is negligible (green and red curves are almost indistinguishable). The perturbation potential due to the $J_3$ term on a $j$ body can be written $U_{J_3} = -\frac{(R)^3 J_3}{r^5} [5 \sin^3(\phi_j) - 3 \sin(\phi_j)]$. Knowing that $\sin(\phi_j) = \sin(i_j) \sin(\nu_j + \omega_j)$ with $i_j, \nu_j$, and $\omega_j$, respectively, the inclination, the true anomaly, and the argument of periapsis of the orbit of $j$. Each term of the potential due to $J_3$ therefore depends on the sine of the inclination, so $i \approx 0 \Rightarrow \nabla U_{J_3} \approx 0$.

The addition of the Sun in the simulations (light blue curve) seems to create perturbations that tend to increase the number of collisions in the short term (approximately equal to 1 yr). The effect of the solar perturbations increases with the semimajor axes of Phobos and Deimos.

<caption>**Table 1**</caption>

Physical Values Used in This Paper

<table><thead><tr><th>Variable</th><th>Value</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td>$J_2$</td><td>$1.96 \times 10^{-3}$</td><td>Ø</td><td>(1)</td></tr><tr><td>$J_3$</td><td>$3.15 \times 10^{-5}$</td><td>Ø</td><td>(1)</td></tr><tr><td>$J_4$</td><td>$-1.54 \times 10^{-5}$</td><td>Ø</td><td>(1)</td></tr><tr><td>$c_{22}$</td><td>$-5.46322410 \times 10^{-5}$</td><td>Ø</td><td>(1)</td></tr><tr><td>$s_{22}$</td><td>$3.15871652 \times 10^{-5}$</td><td>Ø</td><td>(1)</td></tr><tr><td>$\alpha$</td><td>$1.0878850102669404 \times 10^{-7}$</td><td>rad \cdot day^{-1}</td><td>(2)</td></tr><tr><td>$I$</td><td>$25.25797549$</td><td>deg</td><td>(2)</td></tr><tr><td>$\psi$</td><td>$332.6841708$</td><td>deg</td><td>(2)</td></tr><tr><td>$\chi$</td><td>$[0, 2\pi] + 350.89198d$</td><td>deg</td><td>(3)</td></tr></tbody></table>

**Note.** *l*, *ψ*, and *χ* angles are given in B1950. (1) A. Genova et al. (2016); (2) P. Gurfil et al. (2007); (3) M. E. Davies et al. (1983).

<caption>**Table 2**</caption>

Mars Orbital Elements Used for the Simulations with Varying Obliquity

<table><thead><tr><th>Orbital Element</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Semimajor axis</td><td>1.5237</td><td>au</td></tr><tr><td>Eccentricity</td><td>0.09339</td><td>Ø</td></tr><tr><td>Inclination</td><td>0</td><td>rad</td></tr><tr><td>Argument of periapsis</td><td>[0, 2π]</td><td>rad</td></tr><tr><td>Longitude of ascending node</td><td>[0, 2π]</td><td>rad</td></tr><tr><td>True anomaly</td><td>[0, 2π]</td><td>rad</td></tr></tbody></table>

**Note.** The inclination was set to zero, as changing the inclination of the Martian orbit is equivalent to changing its obliquity in a fixed-frame referential centered on Mars. The hooks indicate a range within which a random value is picked.

The precession and nutation of Mars, on the other hand, have no significant effect on collisions. This can be explained by a very slow precession (about $2 \times 10^{-3} \text{deg} \cdot \text{yr}^{-1}$, or 203 deg in 100,000 yr, less than one complete revolution) and a fairly weak nutation on this timescale (oscillation between 23° and 27°).

Finally, taking into account the potential due to the sectorial terms $c_{22}$ and $s_{22}$ has a negligible effect, except for $a_{phobos} = 6R_{Mars}$ (yellow curve). This is the configuration in which Phobos is very close to the synchronous orbit, which most likely induces a resonance that strongly disturbs its orbit. So it seems that the addition of all these physical effects is not sufficient to alter the chances of collisions between Phobos and Deimos after 100,000 yr.

Moreover, as visible in Figure 4, the new configuration does not significantly alter the impact velocity; thus, as in R. Hyodo et al. (2022), the collisions are destructive

## 4.3. The Chaotic Obliquity of Mars

As demonstrated by J. Laskar & P. Robutel (1993), the obliquity of Mars is chaotic and lies between 0 deg and 60 deg. Therefore, it cannot be asserted that when the potential destruction of the progenitor occurred, the rotation axis of Mars was in a configuration similar to the one of 1950. In total, 1000 simulations were therefore carried out for different values of Mars' obliquity, ranging from 0 to 60 deg, by randomly distributing the semimajor axis of Phobos between 5.5 and 6 Mars radii. The orbital parameters of Mars are described in Table 2. As we can no longer use W. R. Ward's (1979) model to compute the revolution of Mars' axis, we now derive it
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Cumulative collision fraction as a function of time for different values of Mars' obliquity. We tried to change Mars' obliquity to see if this could affect the probability of collision of postfragmentation Phobos and Deimos. A total of 1000 simulations were made for each value of the obliquity, and all the physical effects described in this paper were considered. The Martian orbital elements for these simulations are given in Table 2.</figcaption></figure>

<table><caption>**Table 3**
Orbital Elements of the Fragments Used in the Simulations Where More Than Two Fragments Were Considered</caption><thead><tr><td>Orbital Element</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Semimajor axis</td><td>[5, 7.5]</td><td>R<sub>Mars</sub></td></tr><tr><td>Eccentricity</td><td>[0, 0.35]</td><td>Ø</td></tr><tr><td>Inclination</td><td>[0.015, 0.021]</td><td>rad</td></tr><tr><td>Argument of periapsis</td><td>[0, 2π] *</td><td>rad</td></tr><tr><td>Longitude of ascending node</td><td>[0, 2π] *</td><td>rad</td></tr><tr><td>True anomaly</td><td>[0, 2π]</td><td>rad</td></tr></tbody></table>

**Note.** The “*” next to the longitude of ascending node and argument of periapsis values is to notify that, unlike other orbital elements, those values are not independent for each of the fragments. Indeed, we assumed that the fragmentation happened at the node, so the longitude of the ascending node must be the same for all fragments, and the argument of periapsis must be computed in such a way that all orbits cross each other at the node. This last condition is verified when $\omega = \arccos\left(\frac{a(1-e^2)-d_f}{d_f e}\right)$, with a, e, and $\omega$ the semimajor axis, eccentricity, and argument of periapsis of the fragment. $d_f$ is the orbital radius at the fragmentation point, and it is the same for all of the fragments.

directly from the Sun's position and velocity within the simulations.

Looking at Figure 5, we can see that a greater obliquity tends to lower the probability of collision. However, even in the most extreme case, the probability of collision remains greater than 94%. Chaos in Mars' axis of rotation therefore does not alter the results.

### 4.4. Quantity of Fragments

The hypothesis tested by R. Hyodo et al. (2022) is the survival of Phobos and Deimos, assuming they are the only two

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Cumulative fraction of simulations in which all the bodies collided for different numbers of fragments. We made simulations considering all the physical effects presented in this paper, but we considered that the fragmentation led to more than two fragments. The orbital elements repartitions are given in Table 3, and we did 1000 simulations for each number of fragments. We consider that a simulation undergoes “full collision” when all the fragments of the simulation collide with each other.</figcaption></figure>

remnants of the fragmentation. However, if a progenitor did fragment, it is conceivable that there were originally more than two fragments, of which Phobos and Deimos would now be the only survivors. It is of interest to verify if the survival probabilities of two objects are altered if the number of fragments is greater than two. To test this hypothesis, we have generated a number $n$ ($n > 2$) of objects, of which one has the mass of Phobos, a second has the mass of Deimos, and the others have random masses between those of Phobos and Deimos. The orbital elements of all these objects are then randomly distributed within the ranges given by Table 3, which correspond to the maxima and minima of Phobos and Deimos in the results presented by A. Bagheri et al. (2021). To ensure consistency with a fragmentation scenario, all the orbits cross at the ascending node, corresponding to the fragmentation point of the progenitor, which is located between 5.9 and 6.9 Mars radii. No additional constraints are applied to the orbital element. Note that, unlike the two-fragment case, it is not possible to constrain the orbits to have a single point of contact. In order to maximize the fragments' chances of survival, we have considered fusion of fragments after the collision, rather than destructive impacts. The advantage of this assumption is its ability to quickly add up the momentum of the objects, thereby moving them away from each other.

The results of these simulations (Figure 6) show that there is strictly no case in which at least two fragments survive past 100,000 yr.

### 4.5. Tidal Dissipation in Mars, Phobos, and Deimos

To be exhaustive, we also wanted to verify whether the tidal forces could impact the survival rate of the satellites on such short timescales. As migration due to tidal forces is a very slow phenomenon, we did not expect it to have a strong impact on our results. However, the configuration of Mars and its satellites could have been quite different at the time of their formation (hotter Mars, leading to higher $k_2$ value, different values of $Q$ because of the Phobos and Deimos position after
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Cumulative fraction of collisions as a function of simulation time, with different values of tides. We compare the cumulative fraction of collisions in simulations where we included all the effects presented in this paper. To maximize the effect of the tidal forces, we have considered extreme values of $k_2/Q$ significantly above the current values. It appears that even under such extreme circumstances, the collision rate after 25,000 yr remains above 89% in every case.</figcaption></figure>

<table><caption>**Table 4**
Tidal Values Used for the Simulations</caption><thead><tr><td>Body</td><td>k<sub>2</sub></td><td>Q</td></tr></thead><tbody><tr><td>Mars</td><td>0.5</td><td>10</td></tr><tr><td>Phobos</td><td>0.01</td><td>1000</td></tr><tr><td>Deimos</td><td>0.5</td><td>1000</td></tr></tbody></table>

*Note.* The values of $k_2$ and $Q$ are not realistic as each of them were chosen to maximize the effect of the tidal forces as explained in Section 4.5.

their formation, etc.) leading to more important effects of the tidal forces on the satellites orbits than their present values.

To verify if tides could indeed affect the collision rate, we have included them with very high values of $k_2/Q$. We took a value of 0.5 for Mars $k_2$ and 10 for its $Q$, while the current values that are used in A. Bagheri et al. (2021) are, respectively, 0.169 and 95. For the same reasons, we also considered values of $k_2/Q$ for Phobos and Deimos that are 10 times bigger than those used by A. Bagheri et al. (2021). All those values are available in Table 4. As we are assuming to be just after a potential fragmentation, we can not assume synchronicity between the moons' rotation rate and their revolution angular velocity. For this reason, in each simulation, the initial values of the rotation speed of the moons were randomly chosen between 0 and 10 times Mars' rotation rate. The $c_{22}$ used to compute the change in rotation rate for Phobos and Deimos are, respectively, 0.016476 (K. Willner et al. 2014) and 0.04773 (K. Huang et al. 2024). Integrating the rotation of both Phobos and Deimos was quite demanding in terms of computation time. For this reason, we had to stop the simulation after 25,000 yr.

We can see in Figure 7 that the tides on Mars have very little impact on the collision probability. Unexpectedly, the tides inside the bodies seem to have a stronger effect on the collision probability of simulations in which the satellites are further from Mars. This is actually due to a combination of effects: satellites further from Mars in this case are closer to the synchronous orbit, thus are very perturbed by the triaxial shape

of Mars. This lowering in the collision probability for these values of semimajor axis is thus due to a combination of the tidal forces and the resonance with Mars' rotation. Overall, the extreme tidal forces that we included are barely able to lower the collision probability by a couple of percent after 25,000 yr.

## 5. Conclusion and Discussions

As we showed in Section 4, we could not come up with a configuration that lowers the collision probability below 90%. Scenarios in which Phobos and Deimos could have survived for more than 100,000 yr after their formation account for less than 10% in any of the cases, highly conducive to survival that we considered in this study. Thus, the conclusions of R. Hyodo et al. (2022) remain reliable, despite some weaknesses in the dynamic modeling of the problem. We want to draw the reader's attention to the fact that these results do not in themselves invalidate the hypothesis proposed by A. Bagheri et al. (2021). As presented in Section 1, there is currently no formation scenario of Phobos and Deimos able to account for all the observations and modeling results of the two moons. The Martian Moons eXploration mission, which will launch in the coming years, will provide further insights into the composition of Mars' moons (especially Phobos) and will help clarify whether the satellites have an extra-Martian or in situ origin.

## Acknowledgments

The authors thank Michael Efroimsky for suggesting the idea of testing the results of R. Hyodo et al. (2022) by incorporating additional physical effects and for his advice throughout this work. The authors also thank Pini Gurfil for his continuous guidance and advice. We would also like to thank two anonymous reviewers whose comments allowed the authors to greatly improve the manuscript.

## Data Availability

The input parameter of each simulation used in this paper is available on Zenodo at doi:10.5281/zenodo.14006551.
ORCID iDs

Ryan Dahoumane https://orcid.org/0009-0005-8382-8181

References

Aleshkina, E. Y. 2009, SoSyR, 43, 71

Arrighi, P., & Dowek, G. 2015, in Theory and Practice of Natural Computing (TPNC 2015), 9477, ed. A. H. Dediu, L. Magdalena, & C. Martín-Vide (Berlin: Springer), 137

Bagheri, A., Khan, A., Efroimsky, M., Kruglyakov, M., & Giardini, D. 2021, NatAs, 5, 539

Craddock, R. A. 2011, Icar, 211, 1150

Davies, M. E., Abalakin, V. K., Lieske, J. H., et al. 1983, CeMec, 29, 309

Everhart, E. 1985, IAUC, 83, 185

Genova, A., Goossens, S., LeMoine, F., et al. 2016, Icar, 272, 228

Goldreich, P. 1965, AJ, 70, 5

Gurfil, P., Lainey, V., & Efroimsky, M. 2007, CeMDA, 99, 261

Hesselbrock, A. J., & Minton, D. A. 2017, NatGe, 10, 266

Huang, K., Zhang, L., Yang, Y., Ye, M., & Li, Y. 2024, RemS, 16, 1174

Hunten, D. M. 1979, Icar, 37, 113

Hyodo, R., Genda, H., Sekiguchi, R., Madeira, G., & Charnoz, S. 2022, PSJ, 3, 204

Lainey, V. 2002, PhD thesis, Observatoire de Paris, https://www.imcce.fr/content/medias/publications/publications-recherche/theses-habilitations/docs/Lainey_These.pdf

Lainey, V., Dehant, V., & Pätzold, M. 2007, A&A, 465, 1075

Lainey, V., Duriez, L., & Vienne, A. 2004, A&A, 420, 1171

Laskar, J., & Robutel, P. 1993, Natur, 361, 608

Madeira, G., Charnoz, S., Zhang, Y., et al. 2023, AJ, 165, 161

Rein, H., & Spiegel, D. S. 2014, MNRAS, 446, 1424

Rosenblatt, P. 2011, A&ARv, 19, 44

Singer, S. F. 2007, in First Int. Conf. on the Exploration of Phobos and Deimos, 1377 (Houston, TX: Lunar and Planetary Institute), 36

Ward, W. R. 1979, JGR, 84, 237

Willner, K., Shi, X., & Oberst, J. 2014, P&SS, 102, 51

Yoder, C. F. 1995, Astrometric and Geodetic Properties of Earth and the Solar System (Washington, DC: American Geophysical Union)