Numerical model of Phobos’ motion incorporating the effects of free rotation

We developed the equations of motion in a planetocentric (Mars) reference frame with fixed axes align with the ICRF. Hence, in such a reference frame, the position vectors of the Sun and planets relative to Mars can be easily retrieved from the ephemerides. In this paper, we use INPOP19a (Fienga et al., 2019), the newest version of the Intégrateur numérique planétaire de l’Observatoire de Paris (INPOP) planetary ephemerides provided by Institut de Mécanique Céleste et de Calcul des Éphémérides (IMCCE) to get the position and velocity vectors of the Sun and planets relative to Mars in Barycentric coordinates. This ephemeris incorporates updated the Mars observational data by adding the latest data such as ESA’s MEX observations up to 2017, NASA’s MGS, Mars Odessey (MO), and Mars Reconnaissance Orbiter (MRO) observations available from 1999 to 2014 (Fienga et al., 2020).

The orbital motion of Phobos around Mars in our model is described by using position ${\bf r}\equiv(\mathrm{r_{x}},\mathrm{r_{y}},\mathrm{r_{z}})$ and velocity ${\bf v}\equiv(\mathrm{v_{x}},\mathrm{v_{y}},\mathrm{v_{z}})$ in Cartesian coordinates. We can easily write the classical differential equations of relative motion in the planetocentric system as:

where ${\bf F}_{\mathrm{p}}$ and ${\bf F}_{\mathrm{0}}$ indicate the whole external forces exerted on Phobos and Mars, $\mathrm{m_{p}}$ is the mass of Phobos, $\mathrm{m_{0}}$ is the mass of Mars, and $t$ is the time expressed in TDB (Barycentric Dynamical Time) timescale.

The forces that induce the relative motion can be divided into a two-body force and a perturbation force made up of three parts: a part for third-body perturbation, a part for tidal perturbation, and a part for relativistic perturbation. Hence, $\frac{d^{2}{\bf r}}{dt^{2}}$ can be rewritten as :

In the barycentric coordinates system, the two-body force contains interaction between Phobos and Mars considered as point mass and interaction between Phobos and Mars, taking the form:

where $G$ is the gravitational constant, $\mathrm{U_{p}}$ and $\mathrm{U_{0}}$ denote the gravity potential induced on the punctual body Mars and Phobos, and $\mathrm{r}$ denotes the norm of $\bf{r}$. The gravity potential can be conveniently represented as a spherical harmonic expansion with normalized coefficients ($\bar{C}_{lm}$, $\bar{S}_{lm}$), and is given by (Kaula, 1966):

where $l$ is the degree, $m$ is the order, $\bar{P}_{lm}$ are the fully normalized Legendre polynomials, $R$ is the reference radius of central body, and ($r,\phi_{lat},\lambda$) are the radius, latitude, and longitude in body-fixed coordinates of central body (where $\phi_{lat}$ for the latitude to distinguish from Euler angle $\phi$ below), respectively.

The value of ${\bf a}_{\mathrm{two-body}}$ in Eq. 2 does not refer to an inertial or Newtonian coordinate system, but it is itself subject to an acceleration by the third body (Brouwer & Clemence, 2013). Hence, denoting the position vector, $\bf{r_{j}}$, of a third-body, $\mathrm{P_{j}}$, (the Sun, a perturbing planet or Deimos) relative to Mars, and the third-body perturbation can be expressed as

Here $\mathrm{m_{j}}$ is the mass of the third-body $\mathrm{P_{j}}$, and $\mathrm{r_{0j}}=|\bf{r_{j}-r}|$. The modeled dynamics included the third-body perturbations due to Deimos, the Earth, Moon, Jupiter, Saturn, and the Sun.

It was Sharpless (1945) who first found that including an acceleration term can improve his fit to the observed longitudes of Phobos. He concluded that it is a tidal acceleration, but with some uncertainty. Finally, Shor (1975) found an aptly calculated value when he used new observation data in the fitting and eventually confirmed the significant influence of the Martian tidal deformation on Phobos’ orbital motion. In principle, the gravitational attraction of both satellites, the planets and the Sun all raise tides on Mars. The elastic response of Mars to these tides is described by Love numbers. The second degree tides are the strongest, and the change in the gravitational potential due to these tides is given by the potential Love number, $k_{2}$. The tidal potential, $V_{tide}$, caused by a single disturbing body is given by (Goossens & Matsumoto, 2008):

where $\mathrm{m_{i}}$ is the disturbing body, $R_{i}$ is the distance between the disturbing body and Mars, $a_{m}$ is the equatorial radius of Mars, $\hat{\bf R}_{i}$ is the unit vector from the centre of Mars to the disturbing body, and $\hat{\bf r}$ is the unit vector from the centre of Mars to observation point (Phobos in this case).

As for the tide raised by Phobos, instead of Eq. 6, we refer to Lainey et al. (2007) for a full description. Thus, the tidal force ${\bf F}_{T}$ acting on Phobos takes the form:

where $\bf\Omega$ is the Martian angular velocity vector and $\Delta t$ is the time delay due to the inelastic response of Mars (Mignard, 1980; Lainey et al., 2007). By taking Eqs. 6 and  7 into account in the calculation of tidal forces on Phobos, the value of ${\bf a}_{\mathrm{tide}}$ can be modeled.

Finally, we introduce the relativistic effects, ${\bf a}_{\mathrm{rel}}$. This model is originally inspired by Einstein-Infeld-Hoffmann (EIH) equations of motion used in the planetary and lunar ephemerides (Folkner et al., 2014; Pitjeva & Pavlov, 2017; Viswanathan et al., 2017). For each body (either Mars or Phobos), $A$, the relativistic acceleration due to interaction with other point masses in solar system, ${\bf a}_{A,rel}$ is given by (Folkner et al., 2014):

where $\beta$ and $\gamma$ are the parametrized post-Newtonian (PPN) parameters and both fixed at 1, with $c$ being the speed of light in vacuum.

Carrying out these formulas Eqs. 3 - 8 for ${\bf a}_{\mathrm{two-body}}$, ${\bf a}_{\mathrm{third-body}}$, ${\bf a}_{\mathrm{tide}}$, and ${\bf a}_{\mathrm{rel}}$ in turn and adding them together, we obtained the orbital motion equations in the planetocentric coordinates.

In this paper, Phobos is modeled as an elastic body according to our previous work on Phobos’ libration (Yang et al., 2020). The orientation of Phobos is integrated from the differential equations for its angular velocities. The angular momentum vectors of Phobos are the product of the angular velocities and the moments of inertia. The angular momentum vectors change with time due to torques and distortion of the body.

To express the variations of Phobos’ orientation in the inertial frame, it would be convenient to define Phobos’ frame aligned with the principal axes of Phobos. Then the orientation of Phobos’ frame with respect to the inertial frame is determined by three Euler angles: $\phi$, $\theta$, and $\psi,$ which evolve over time. The transformation from body-fixed frame (principal axes) to the inertial frame is given by the matrix:

where the rotation matrices, $R_{x}$ and $R_{z}$, are right-handed rotations around the x-axis and z-axis, respectively. Hereinafter, the argument time, $t,$ will be omitted for simplicity.

The instantaneous rate of Euler angle at time, $t$, $\frac{d\phi}{dt}$, $\frac{d\theta}{dt}$, and $\frac{d\psi}{dt}$, along with the Euler angles, can be used to describe ${\bm{\omega}}(t)\equiv(\omega_{1},\omega_{2},\omega_{3})^{\mathrm{T}}$, the angular velocity of Phobos (Goldstein, 2011):

where $\mathrm{T}$ means transpose matrix. Then differentiating Eq. 10 with respect to time we obtain a system of equations linear in $\frac{d^{2}\phi}{dt^{2}}$, $\frac{d^{2}\theta}{dt^{2}}$, and $\frac{d^{2}\psi}{dt^{2}}$, for which we algebraically solve, obtaining:

In this paper, the first- and second- order time derivatives of $X$ are also denoted by $\dot{X}$ and $\ddot{X}$, respectively.

In a rotating system, the change in angular velocity $\bm{\omega}$ is related to torques $\bf N$ and governed by Euler’s rotation equations:

where $\mathbf{I}$ is the moment of inertia tensor. This leads to an equation for $\bm{\omega}$, therefore,

Williams et al. (2001) first gave the numerical formulas for calculating the inertia tensor of a celestial body with elastic property, these formulas were then implemented to generate the lunar ephemerides (Folkner et al., 2014; Pitjeva & Pavlov, 2017) and numerically study the rotation of Phobos (Yang et al., 2020). In this paper, we refer to Yang et al. (2020) for full description to construct the inertia tensor of Phobos. The inertia tensor of the Phobos mainly contains a rigid part $I_{rigid}$, while it is also subject to tidal distortion from Mars and spin distortion via:

Here, $I_{tide}$ and $I_{spin}$ are the inertia tensor due to tidal deformation owing to the attraction of Mars and spin distortion, respectively.

The torques, $\bf N,$ referred to the Phobos’ frame are generally split into two parts: a torque from point-mass (body) A to the Phobos’ figure and a torque from the martian oblateness to the Phobos’ figure:

The detailed description of $\bf I$ and $\bf N$ can refer to Williams et al. (2001), Folkner et al. (2014), Rambaux et al. (2012) and Yang et al. (2020).

The instantaneous state of Phobos’ orientation can be defined completely by six quantities if we choose the Euler angles $(\phi,\theta,\psi)$ and their rates $(\frac{d\phi}{dt},\frac{d\theta}{dt},\frac{d\psi}{dt})$ to describe the state of Phobos’ orientation, Eqs. 10, 11, and 13 form the differential equations for Phobos’ rotational motion.

The orbital and rotational equations of motion developed in the previous two sections are non-linear and cross-coupled, and they are not easily solved to determine the position and Euler angles as function of time. Thus, in the implementation of our dynamical model of Phobos’ motion, we numerically integrated them simultaneously via Runge-Kutta-Fehlberg (RKF) adaptive procedure method (Simos, 1993; Yang et al., 2018). Numerical experiments show that the computational procedure is of good accuracy by taking a fixed step size (10 minutes) during the integration.

The RKF integrator is configured for approximating the solution of the first-order differential equation $y^{\prime}=f(x,y)$ with an initial condition $y(x_{0})=c$, so we rewrite the three second-order orbital equations of motion (Eq. 1) as a set of six first-order differential equations. The result is:

Similarly, the three second-order rotational equations of motion can be rewritten as:

where $\Phi(t)\equiv(\phi,\theta,\psi)^{\mathrm{T}}$.

The six orbital equations of motion are integrated in parallel with six rotational equations of motion (Eqs. 16 and 17). The initial condition state vectors are the three Euler angles and their rates, and position and velocity of Phobos relative to Mars at an initial epoch (J2000.0 in this work). In addition, to control the efficiency and of the numerical integrator and also to validate the computation of the equations of motion and rotation, we use the conservation of energy of the system, that is, neglecte the Martian tidal dissipation. The full list of physical models and parameters used in the dynamical model is given in Table. 1.

There are a number of parameters in our numerical model whose values affect the orbit significantly, such as the Martian gravity field, Phobos’ initial position and velocity, and spin libration of Phobos (Lainey et al., 2007). To reveal the difference between fully coupled approach with the simpler one used so far, we first quote the model now employed in the ephemerides (Jacobson, 2010; Jacobson & Lainey, 2014; Lainey et al., 2020b), namely: the simple model. Since Phobos’ rotation period match its orbit period and in synchronous rotation, assuming that the satellite’s pole normal to its orbit plane, the angle between the satellite’s axis of minimum principal moment of inertia and direction from satellite to Mars is small. In this case, it can be approximated as $\theta=(2e+\mathcal{A})\sin M$, where $e$, $\mathcal{A}$ and $M$ are the satellite’s orbital eccentricity, the libration amplitude (Rubincam et al., 1995), and the mean anomaly of satellite in its orbit, respectively. The force on Mars exerted by Phobos based on this assumption can be calculated as:

where ${\mu_{p}}$ is the $GM$ of Phobos, ${a_{p}}$ is the equatorial radius of Phobos, $C_{20}$ and $C_{22}$ are the second degree harmonic of Phobos, $\hat{\bf{r}}$ is the unit vector directed from Mars toward satellite, and $\hat{\bf{t}}$ is the unit vector in the satellite’s orbit plane normal to $\hat{\bf{r}}$ and in the direction of its orbital motion. Easily, the reactive force acting on the satellite is:

with $\mu_{m}$ denoting the $GM$ of Mars. Using the model presented above, we simulated the current simpler model, but used our own selected physical parameters in Table. 1. Then, we used the six initial conditions (position and velocity) as our adjustment parameters to fit the current ephemeris of Phobos.

We chose the most recent Mars moon ephemerides NOE-4-2020 (Lainey et al., 2020b) as the ”observations” and fit our model to them. The adjustment was carried out in Cartesian planetocentric coordinates J2000 using a sample of 3650 points with a step size of one day (ten years). We integrated over ten years forth starting at the initial epoch Julian day 2451545.0 (J2000.0, TDB timescale). Figure 1 shows the residuals after applying the least-square fitting procedure. The deviations in the position are probably explained by the different physical parameters (such as the Martian dissipation factor, Q, the physical libration, $\mathcal{A}$, the $C_{20}$ and $C_{22}$ of Phobos, etc.) in these two models. This fitting result gives us an optimal reference for studying the differences between the newly full model and the simpler one used so far.

To explore the difference between our full 3D dynamical model and the previous simple model (Jacobson, 2010; Jacobson & Lainey, 2014; Lainey et al., 2020b). We adjusted the full rotational model to the reference ephemeris integrated in Section  4.1. Firstly, the gravity coefficients of Phobos are truncated to the second degree, so as to distinguish the differences that are only due to the modeling.

In this paper, we selected the twelve initial conditions (position, velocity, Euler angles, and their rates) as our adjustment parameters to fit the current reference results integrated from the simulated simple model in previous section. The initial position and velocity are retrieved from NOE-4-2020 directly. Phobos’ Euler angles and its rates refer to Archinal et al. (2018) and the fitting strategy is the same as the previous one. Figures 2 and  3 show the residuals and difference after applying the least-square fitting procedure to our simulated simple model, the deviations of the position are probably explained by the newly introduced latitudinal libration and different longitudinal libration of Phobos in the full model. The evolution of the Euler angles defined respect to the inertial reference frame system for 10 years from 1 Jan 2000 to 2010 are plotted in Fig. 4.

To analyze the rotation behavior, Fig. 5 presents the angles between the direction from Phobos to Mars and Phobos’s axis of minimum principal moment of inertia in the two models, respectively. The differences between the two dynamical models are very small, indicating that the simple model can describe the deviation in the longitude direction very well. The Phobos obliquity is shown in Fig. 6, this quantity is neglected in the simple model, so this may be a reason for the different orbits, even though it is an order of magnitude smaller than the angle in the longitude direction. The difference between the two models are then the instantaneous spin pole on the space. On average, the two poles are moving on the same rate that is shown and confirmed in Fig. 7. However, the thickness of the blue points in Fig. 7 indicate the full model has a large oscillation amplitude of around one degree, compared to the simple model which assumes the Phobos pole is aligned with its orbit normal.

The other main purpose of this work is to examine the effect of introducing higher order gravity field coefficients on rotational and orbital motions. However, because the gravity field of Phobos is hard to determine directly from the current data, especially as the coefficients higher than two degrees are now all derived from the shape of Phobos with a homogeneous density hypothesis. Table. 2 provides the numerical values of these coefficients up to the fourth degree and order from a forward modeling method(Yang et al., 2020).

Figure 8 displays the difference after an adjustment about 10 years of our full model with Phobos’ gravity truncated at a degree and order of three to the simulated simple model in Sect. 4.1. One can wonder why the difference in distance after fitting the model that takes into account the third-degree gravity field of Phobos to simulated simple model are more significant than that of a model that takes into account only the second degree. During our adjustment we found $S_{33}$ to be the cause of large orbit difference, that is, if we neglect the $S_{33}$ and the adjustment will be very close to the case where only the second-degree gravity field is considered. Fig. 9 presents the difference in distance after fitting the full model with the third-degree gravity field, but not $S_{33}$ to the simulated simple model of Phobos. The resulting difference is very close to the full model of Phobos with gravity coefficients truncated at the second degree. This result may come from the newly introduced librations by the third-degree harmonic (Borderies & Yoder, 1990; Rambaux et al., 2012), especially with respect to the relatively large value of $S_{33}$, and the tidal orbital expansion may also bring about this signal (Lainey et al., 2020a). In conclusion, this issue needs to be carefully re-studied in the future based on more reliable gravity field coefficients of Phobos.

We compared the results of fitting the full model considering higher degree $(i=3,\cdots,8)$ gravity field coefficients to the simulated simple model. The difference of the post-fitted orbits containing third- and fourth-order gravity field coefficients is plotted in Fig. 10, it indicates that the difference between the two simulations is at the ”meter” level. In addition, we find that the maximum orbit difference does not exceed one meter   even when directly integrating a model containing higher-order gravity field coefficients using the post-fitted initials of a model containing fourth-order gravity field coefficients, as shown in Fig. 11. This suggests that the degree of gravity field coefficients higher than three can be neglected in the full model concerning the current accuracy. However, this conclusion is strongly dependent on observational data and will need to be revisited in the future if lander tracking data on Phobos’ surface or high-precision rotation data become available.

Some uncertainty about the elastic property of Phobos remains. In this work, Phobos was treated as an elastic body where the tidal Love number, $k_{2}$, was referred to by Le Maistre et al. (2013) to construct the numerical model of motion. The value of $k_{2}$ affects the rotation of Phobos (Williams et al., 2001; Yang et al., 2020) and, thus, the orbital motion (Lainey et al., 2007; Jacobson, 2010; Jacobson & Lainey, 2014). Consequently, in order to identify the effect of different $k_{2}$ on Phobos’ orbital motion, we chose various $k_{2}$ to simulate Phobos’ motion.

By including a range of plausible value for $k_{2}\in[1.83\times 10^{-7},5.3\times 10^{-4}]$ (Le Maistre et al., 2013), our new model uses a second-degree gravity field of Phobos and a different $k_{2}$ value fitted to the simulated simple model that without considering $k_{2}$ and then the differences (with and without) are plotted in Fig. 12. As shown in this figure, the differences in orbit are less than $100\leavevmode\nobreak\ m$ in ten years’ integration when $k_{2}\leq 1.0\times 10^{-4}$. This result suggests that it is difficult to effectively determine the $k_{2}$ value of Phobos from the current data.