© The Authors 2024

1 Introduction

The orbital and rotational motions of Phobos have always been of particular interest the keys to understanding the origin and formation of the Martian moons overall (Rosenblatt 2011). Since its discovery in 1877 and thanks to the development of observation technologies and the first Earth-based and subsequent spacecraft observations, a range of ephemerides, including analytical theories and numerical theories, have been developed (Jacobson & Lainey 2014).

One of the first analytical theories for the description of the satellites’ orbital motion around Mars was developed by Struve (1911). In this theory, Phobos and Deimos are restricted to an elliptical orbit. To explain the satellites’ secular perturbation due to the Mars oblateness and the Sun, Struve introduced the precession of pericenter longitude, ϖ, and ascending node, Ω.

However, since the perturbation of the orbits are conics, the satellites’ semi-major axes and mean motions do not satisfy Kepler’s third law very well (Jacobson & Lainey 2014). Shor (1975) revised this theory through adding an acceleration term to the mean orbital longitudes to explain the tidal acceleration and estimated the parameters based on observational data collected from 1877 to 1973. Using the variation of arbitrary constants method (Brouwer & Clemence 2013), Sinclair (1978, 1989) constructed a dynamical theory to obtain an approximate analytical solution to the satellites’ equations of orbital motion. The elements (a, e, λ, ϖ, Ω) and the sine of the inclination of the satellite’s orbit to the Laplace plane are corrected with secular and periodic terms. Besides the elements, the mean longitude accelerations, the mass of Mars, and the zonal harmonics of the Martian gravity field up to the fourth degree were also included in this theory. To improve the accuracy of the orbit of Phobos and support Soviet Phobos 2 mission, Morley (1989) optimized Sinclair’s theory by taking into account more terms including some of second-order perturbations.

Given that the Martian satellites are small, it was found that theories of artificial satellites of the Earth could be adapted as models for their orbits. Unlike the theories of Struve or Sinclair, those based on the artificial satellite method used the Mars equator, but not the Laplace plane, as the reference plane. With the advent of artificial satellite theory and digital computer technology, a number of analytical and semi-analytical theories have been developed to describe their orbital motion (Ivanov et al. 1988; Kudryavtsev 1993; Emelyanov & Nasonova 1989; Emelyanov et al. 1993).

The most well known and also high-precision semi-analytical approach is Chapront-Touzé theory (Chapront-Touzé 1988, 1990a,b). It takes into account the perturbative effects of the Sun and planets, the Martian gravity field, the mutual perturbations of Phobos and Deimos, the precession and nutation of the Martian equator, and the perturbations due to the gravity field of Phobos. Unlike the orbital elements used in previous theory, the satellite Cartesian coordinates and their velocities in the Mars equator of date coordinate system is first adapted into the model. The most significant improvement of this theory is that it models the libration of Phobos based on the work of Borderies & Yoder (1990). However, due to the semi-analytical nature of the theory, it is difficult to revise according to new observations or new values for fundamental constants such as the coefficients of Martian and Phobos’ gravity field, that is, the numerical terms would have to be regenerated. There have been no updates to the theory published since.

Another approach to studying the motion of Phobos is based on a purely numerical method. The first elegant, purely numerically integrated ephemeris for the Martian satellites was constructed by Lainey et al. (2007) based on their software called Numerical Orbit and Ephemerides (NOE; Lainey et al. 2004a,b). They modeled the perturbative effects of the Sun, Jupiter, Saturn, Earth, and Moon using DE406 ephemeris, an aspherical Martian gravity field up to the tenth degree and order (Tyler et al. 2003), mutual perturbations among the two satellites, effects of tides raised on Mars by Phobos (Mignard 1980), IAU2000 Martian precession and rotation (Seidelmann et al. 2002), as well as the C₂₀ and C₂₂ of Phobos. Lainey’s integration was implemented in planetocentric Cartesian coordinates referred to the International Celestial Reference Frame (ICRF) with an integrator called RA15 (Everhart 1985). They also fitted to observations from 1877 to 2005 and included the spacecraft observations of Mars Global Surveyor (MGS), and Mars EXpress (MEX). The ultimately estimated position accuracy was roughly one kilometer or more with respect to Mars over a century (Lainey et al. 2007). To support the planned Russian Phobos-Grunt mission (Marov et al. 2004), Shishov (2008) employed a similar force model but using Lagrange’s equations to integrate the orbital elements. A numerical integration method of the eighth order of accuracy (Stepaniants & Lvov 2000) was used.

In Lainey’s model, they realized that the spin librations of Phobos would greatly affect the accuracy of the model. While the Moon revolves around the Earth in a synchronous rotation state, Phobos is in a synchronous spin-orbit resonance around Mars. Thus, these authors forced Phobos in a 1:1 spin-orbit resonance directly in their model. Jacobson (2010) further developed this model by assuming that Phobos is in synchronous rotation and taking into account only the longitudinal libration (amplitude about 1.24° by Willner et al. 2014) during its orbital motion around Mars. It was perfectly reasonable to make such an assumption at the time, since there were no observations of Phobos’ Euler angles can be used easily.

Fortunately, this situation will be changed by future missions, which include Martian Moons eXploration (MMX) to be launched in 2024 (Kawakatsu et al. 2017; Usui et al. 2018). A spacecraft dedicated to the satellite will give us a clearer understanding of the shape and physical characteristics such as gravity field and rotational motion of Phobos in the future.

To prepare for such future missions, we developed a new pure numerical model of Phobos’ motion that takes no assumption of Phobos’ rotation. Phobos’ rotational motion is firstly introduced by Euler’s rotational equations (Cappallo et al. 1981; Rambaux et al. 2012; Yang et al. 2020). This full 3D model will enable us to simulate the motion of Phobos with a more realistic case based on the future observation data. To address this issue, we present a revised dynamical model of Phobos’ motion around Mars in this paper. Section 2 describes the equations of Phobos’ orbital and rotational motion in an appropriate reference system. Section 3 presents the variational equations of the rotational motion. In Sect. 4, we use adjustment test simulations to compare the difference between our model and the simple model currently in use. Finally, we briefly discuss the results and present the main conclusions of this paper in Sect. 5.

2 Dynamical model

The first step to numerically integrating the motion of Phobos is introducing the modeling process. We naturally divide the modeling process into orbit modeling and rotation modeling.

2.1 Translational equations of motion

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 Celeste 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 r ≡ (r_x, r_y, r_z) and velocity v ≡ (v_x, v_y, v_z) in Cartesian coordinates. We can easily write the classical differential equations of relative motion in the planetocentric system as: $$\frac{{\text{d}}^{2}r}{\text{d}{t}^2}=\frac{F_{\text{p}}}{m_{\text{p}}}-\frac{F_0}{m_0},$$(1)

where F_p and F₀ indicate the whole external forces exerted on Phobos and Mars, m_p is the mass of Phobos, m₀ 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{{\text{d}}^{2}r}{\text{d}{t}^2}$ can be rewritten as: $$\frac{{\text{d}}^{2}r}{\text{d}{t}^2}=a_{\text{two-body }}+a_{\text{third-body }}+a_{\text{tide }}+a_{rel}\text{. }$$(2)

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: $$a_{\text{two–body }}=-\frac{G\left(m_{\text{p}}+m_0\right)r}{r^3}-G{m}_0\nabla U_{\text{p}}+G{m}_0\nabla {\text{U}}_0,$$(3)

where G is the gravitational constant, U_p and U₀ denote the gravity potential induced on the punctual body Mars and Phobos, and r denotes the norm of r. The gravity potential can be conveniently represented as a spherical harmonic expansion with normalized coefficients $\left({\overline{C}}_{lm},{\overline{S}}_{lm}\right)$, and is given by (Kaula 2013): $$\begin{array}{l}U\left(r,{\varphi }_{\text{lat }},\lambda \right)=\frac{1}{r}\times {\displaystyle \sum _{l=0}^{l_{\max }}}{\left(\frac{R}{r}\right)}^l\hfill \\ \times {\displaystyle \sum _{m=0}^l}{\overline{P}}_{lm}\left(\sin {\varphi }_{\text{lat}}\right)\left[{\overline{C}}_{lm}\cos m\lambda +{\overline{S}}_{lm}\sin m\lambda \right],\hfill \end{array}$$(4)

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

The value of a_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, r_j, of a third-body, P_j, (the Sun, a perturbing planet or Deimos) relative to Mars, and the third-body perturbation can be expressed as $$a_{\text{third–body }}={\displaystyle \underset{j}{{\sum }^{\text{}}}}G{m}_j\left(\frac{r_j-r}{r_{0j}^3}-\frac{r_j}{r_j^3}\right).$$(5)

Here m_j is the mass of the third-body P_j, and r_0j·= |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₂. The tidal potential, V_tide, caused by a single disturbing body is given by (Goossens & Matsumoto 2008): $$V_{\text{tide }}=\frac{k_2}{2}\frac{G{m}_i}{R_i^3}\frac{a_m^5}{r^3}\left(3{\left({\widehat{R}}_i\cdot \widehat{r}\right)}^2-1\right),$$(6)

where 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, ${\widehat{R}}_i$ is the unit vector from the centre of Mars to the disturbing body, and $\widehat{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 F_T acting on Phobos takes the form: $$F_T=-\frac{3k_{2}G{m}_{\text{p}}^2{a}_m^5}{r^8}\left(r+\text{Δ}t\left[\frac{2r(r\cdot v)}{r^2}+(r\times \Omega +v)\right]\right),$$(7)

where Ω is the Martian angular velocity vector and ∆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 a_tide can be modeled.

Finally, we introduce the relativistic effects, a_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, a_A,rel is given by (Folkner et al. 2014): $$\begin{array}{l}\begin{array}{c}{a}_{A,\text{ rel }}={\displaystyle \underset{j\ne i}{{\sum }^{\text{}}}}\frac{G{m}_j\left(r_j-r_i\right)}{r_{ij}^3}\{-\frac{2(\beta +\gamma )}{c^2}{\displaystyle \underset{k\ne i}{{\sum }^{\text{}}}}\frac{G{m}_k}{r_{ik}}-\frac{2\beta -1}{c^2}\frac{G{m}_k}{r_{jk}}\\ +\gamma {\left(\frac{v_i}{c}\right)}^2+(1+\gamma ){\left(\frac{v_j}{c}\right)}^2-\frac{2(1+\gamma )}{c^2}{v}_i\cdot v_j\end{array}\\ \begin{array}{l}-\frac{3}{2c^2}{\left[\frac{\left(r_i-r_j\right)\cdot v_j}{r_{ij}}\right]}^2+\frac{1}{2c^2}\left(r_i-r_j\right)\cdot \frac{\text{d}{v}_j}{\text{d}t}\}\hfill \\ +\frac{1}{c^2}{\displaystyle \underset{j\ne i}{{\sum }^{\text{}}}}\frac{G{m}_j}{r_{ij}^3}\left\{\left(r_i-r_j\right)\cdot \left[(2+2\gamma )v_i-(1+2\gamma )v_j\right]\right\}\left(v_i-v_j\right)\hfill \\ +\frac{3+4\gamma }{2c^2}{\displaystyle \underset{j\ne i}{{\sum }^{\text{}}}}\frac{G{m}_j}{r_{ij}}\frac{\text{d}{v}_j}{\text{d}t},\hfill \end{array}\end{array}$$(8)

where β and γ 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 a_two–bodý, a_third–body, a_tide, and a_rei in turn and adding them together, we obtained the orbital motion equations in the planetocentric coordinates.

2.2 Evolution of Phobos’ orientation

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: ϕ, θ, and ψ, which evolve over time. The transformation from body-fixed frame (principal axes) to the inertial frame is given by the matrix: $$r_{\text{ICRF}}=R_z(-\varphi (t))R_x(-\theta (t))R_z(-\psi (t))r_{\text{PA}},$$(9)

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

The instantaneous rate of Euler angle at time, t, $\frac{\text{d}\varphi }{\text{d}t},\frac{\text{d}\theta }{\text{d}t}$, and $\frac{\text{d}\psi }{\text{d}t}$, along with the Euler angles, can be used to describe ω(t) ≡ (ω₁, ω₂, ω₃)^T, the angular velocity of Phobos (Goldstein 2011): $$\omega =\left[\begin{array}{ccc}\sin \theta \sin \psi & \cos \psi & 0\\ \sin \theta \cos \psi & -\sin \psi & 0\\ \cos \theta & 0& 1\end{array}\right]\left[\begin{array}{c}\frac{\text{d}\varphi }{\text{d}t}\\ \frac{\text{d}\theta }{\text{d}t}\\ \frac{\text{d}\psi }{\text{d}t}\end{array}\right],$$(10)

where T means transpose matrix. Then differentiating Eq. (10) with respect to time we obtain a system of equations linear in $\frac{{\text{d}}^2\varphi }{\text{d}{t}^2},\frac{{\text{d}}^2\theta }{\text{d}{t}^2}$, and $\frac{{\text{d}}^2\psi }{\text{d}{t}^2}$, for which we algebraically solve, obtaining: $$\left[\begin{array}{c}\frac{{\text{d}}^2\varphi }{\text{d}{2}^2}\\ \frac{{\text{d}}^2\theta }{\text{d}{2}^2}\\ \frac{{\text{d}}^2\psi }{\text{d}{t}^2}\end{array}\right]=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\text{d}\omega }{\text{d}t}+\left[\begin{array}{c}\frac{\dot{\theta }(\dot{\psi }-\dot{\varphi }\cos \theta )}{\sin \theta }\\ -\dot{\varphi }\dot{\psi }\sin \theta \\ -\ddot{\varphi }\cos \theta +\dot{\varphi }\dot{\theta }\sin \theta \end{array}\right].$$(11)

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 ω is related to torques N and governed by Euler’s rotation equations: $$\frac{\text{d}}{\text{d}t}(I\omega )+\omega \times I\omega =N,$$(12)

where I is the moment of inertia tensor. This leads to an equation for ω, therefore, $$\frac{\text{d}\omega }{\text{d}t}=I^{-1}\left(N-\frac{\text{d}I}{\text{d}t}\omega -\omega \times I\omega \right).$$(13)

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: $$I=I_{\text{rigid }}+I_{\text{tide }}+I_{\text{spin }}\text{. }$$(14)

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, 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: $$N=N_{\text{Pho,fig-pm }}+N_{\text{Pho,fig-fig }},$$(15)

The detailed description of I and 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 (ϕ, θ, ψ) and their rates $\left(\frac{\text{d}\varphi }{\text{d}t},\frac{\text{d}\theta }{\text{d}t},\frac{\text{d}\psi }{\text{d}t}\right)$ to describe the state of Phobos’ orientation, Eqs. (10), (11), and (13) form the differential equations for Phobos’ rotational motion.

2.3 Solution of the equations by numerical integration

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 min) during the integration.

The RKF integrator is configured for approximating the solution of the first-order differential equation y′ = f (x, y) with an initial condition y(x₀) = 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: $$\begin{array}{c}\frac{\text{d}r}{\text{d}t}=v\\ \frac{\text{d}v}{\text{d}t}=\frac{{\text{d}}^{2}r}{\text{d}{t}^2}.\end{array}$$(16)

Similarly, the three second-order rotational equations of motion can be rewritten as: $$\begin{array}{c}\frac{\text{dΦ}}{\text{d}t}=\dot{\text{Φ}}\\ \frac{\text{d}\dot{\text{Φ}}}{\text{d}t}=\frac{{\text{d}}^2\text{Φ}}{\text{d}{t}^2},\end{array}$$(17)

where Φ(t) ≡ (ϕ, θ, ψ)^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, neglected the Martian tidal dissipation. The full list of physical models and parameters used in the dynamical model is given in Table. 1.

3 Partial derivatives

In order to combine our new dynamical model with future higher-precision observations (e.g. MMX) to further the study of Phobos, we rely upon the iterative use of a linear least-squares estimator. Thus, in addition to the model for motion, it is necessary to generate partial derivatives of the state vectors with respect to the rotational and orbital motion initial conditions and parameters we are interested in at all times.

The partial derivatives of the state vectors with respect to initial position and velocity of Phobos and dynamical parameters have been thoroughly studied (Peters 1981; Taylor 1998; Montenbruck et al. 2002; Lainey et al. 2004a,b). Here, we focus on the six initial conditions of the rotational state vectors. Denoting the state vectors as $X(t)={(r(t),\stackrel{}{r\dot{}}(t),\text{Φ}(t),\dot{\text{Φ}}(t))}^T$, the state transition matrix could then be obtained from: $$\frac{\text{d}}{\text{d}t}\left(\frac{\partial X(t)}{\partial X\left(t_0\right)}\right)=\frac{\partial \dot{X}(t)}{\partial X(t)}\left(\frac{\partial X(t)}{\partial X\left(t_0\right)}\right),$$(18)

and the initial value, $${\left(\frac{\partial X(t)}{\partial X\left(t_0\right)}\right)|}_{t=t_0}=1_{12\times 12}.$$(19)

We can write $\frac{\partial \dot{X}(t)}{\partial X(t)}$ explicitly, $$\begin{array}{l}\frac{\partial \stackrel{.}{X}(t)}{\partial X(t)}=\left(\begin{array}{llll}\partial \stackrel{.}{r}/\partial r\hfill & \partial \stackrel{.}{r}/\partial \stackrel{.}{r}\hfill & \partial \stackrel{.}{r}/\partial \text{Φ}\hfill & \partial \stackrel{.}{r}/\partial \dot{\text{Φ}}\hfill \\ \partial \stackrel{\mathrm{..}}{r}/\partial r\hfill & \partial \stackrel{\mathrm{..}}{r}/\partial r\hfill & \partial \stackrel{\mathrm{..}}{r}/\partial \text{Φ}\hfill & \partial \stackrel{\mathrm{..}}{r}/\partial \dot{\text{Φ}}\hfill \\ \partial \dot{\text{Φ}}/\partial r\hfill & \partial \dot{\text{Φ}}/\partial \stackrel{.}{r}\hfill & \partial \dot{\text{Φ}}/\partial \text{Φ}\hfill & \partial \dot{\text{Φ}}/\partial \dot{\text{Φ}}\hfill \\ \partial \ddot{\text{Φ}}/\partial r\hfill & \partial \ddot{\text{Φ}}/\partial \stackrel{.}{r}\hfill & \partial \ddot{\text{Φ}}/\partial \text{Φ}\hfill & \partial \ddot{\text{Φ}}/\partial \text{Φ}\hfill \end{array}\right)\\ =\left(\begin{array}{cccc}0& 1& 0& 0\\ \partial \stackrel{\mathrm{..}}{r}/\partial r& \partial \stackrel{\mathrm{..}}{r}/\partial r& \partial \stackrel{\mathrm{..}}{r}/\partial \text{Φ}& \partial \stackrel{\mathrm{..}}{r}/\partial \text{Φ}\\ 0& 0& 0& 1\\ \partial \ddot{\text{Φ}}/\partial r& \partial \ddot{\text{Φ}}/\partial \stackrel{.}{r}& \partial \ddot{\text{Φ}}/\partial \text{Φ}& \partial \ddot{\text{Φ}}/\partial \text{Φ}\end{array}\right).\end{array}$$(20)

To obtain the expression of none-zero components in Eq. (20), we split them up into four cases:

• 1.

  $\partial \stackrel{}{r¨}/\partial r$ and $\partial \ddot{r}/\partial \dot{r}$ have been thoroughly studied (Peters 1981; Taylor 1998; Montenbruck et al. 2002; Lainey et al. 2004a,b) and were thus omitted here. However, it should be noticed that $\partial \ddot{r}/\partial r$ and $\partial \ddot{r}/\partial \dot{r}$ are a little bit different with previous results since the rotational motion of Phobos were described by Euler angles here. By using Eq. (9), we get: $$\xi =R_z(\psi )R_x(\theta )R_z(\varphi )r\equiv ℛ(\text{Φ})r,$$(21)

  where ξ denotes the coordinates referred to the body-fixed (principal axis) system, if we differentiate Eq. (21) with respect to time, we obtain the relationship of the velocity in two systems: $$\frac{\text{d}\xi }{\text{d}t}=\frac{\partial ℛ(\text{Φ})}{\partial t}r+ℛ(\text{Φ})\dot{r}.$$(22)

  Thus, $$\frac{\partial \xi }{\partial r}=ℛ(\text{Φ})=\frac{\partial \dot{\xi }}{\partial \dot{r}}.$$(23)

• 2.

  Based on the Eqs. (2) and (3), we are able to conclude that $\partial \stackrel{}{r¨}/\partial \dot{\text{Φ}}=0$ and $\partial \stackrel{}{r¨}/\partial \text{Φ}=G{m}_0\frac{\partial }{\partial \text{Φ}}\nabla U_{\text{p}}$, respectively. To calculate the $\frac{\partial }{\partial \text{Φ}}\nabla U_{\text{p}}$, the inverse transformation from body-fixed frame (principal axes) of Phobos to the inertial frame is needed. This is discussed in the next section.

• 3.

  To compute $\partial \ddot{\text{Φ}}/\partial r$ and $\partial \ddot{\text{Φ}}/\partial \dot{r}$, we differentiate $\ddot{\text{Φ}}$ with respect to the position and velocity of Phobos relative to Mars in inertial frame. Based on Eqs. (11) and (13) in this paper, it shows that we only need to evaluate $\partial \dot{\omega }/\partial r$ and $\partial \dot{\omega }/\partial \dot{r}$. This problem is also can be divided into two steps, a transformation from body-fixed frame (principal axes) of Phobos to the inertial frame via the Euler angles and a derivates with respect to the position or velocity, the latter was studied in previous literature.

• 4.

  $\partial \ddot{\text{Φ}}/\partial \text{Φ}$ and $\partial \ddot{\text{Φ}}/\partial \dot{\text{Φ}}$ will be described in detail in this work.

With the above analysis, we can conclude that the key to establishing the variational equation are $\partial \ddot{\text{Φ}}/\partial \text{Φ,\hspace{0.17em}}\partial \ddot{\text{Φ}}/\partial \dot{\text{Φ}}$, and transformation matrix between body-fixed frame (principal axes) of Phobos and the inertial frame. This paper gives these differential equations explicitly.

Referring back to the defining Eq. (11), we differentiate it for the φ with respect to each component of the state vectors: $$\frac{\partial \ddot{\text{Φ}}}{\partial \varphi }=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \varphi }+\left[\begin{array}{c}0\\ 0\\ \frac{\partial \ddot{\varphi }}{\partial \varphi }\cos \theta \end{array}\right],$$(24) $$\frac{\partial \ddot{\text{Φ}}}{\partial \theta }=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \theta }+\left[\begin{array}{c}(\dot{\varphi }\dot{\theta }-\ddot{\varphi })\cot \theta \\ -\varphi \dot{\psi }\cos \theta \\ \ddot{\varphi }\sin \theta +\cos \theta \left(\dot{\varphi }\dot{\theta }-\frac{\partial \ddot{\varphi }}{\partial \theta }\right)\end{array}\right]\text{, }$$(25) $$\frac{\partial \ddot{\text{Φ}}}{\partial \psi }=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \psi }+\left[\begin{array}{c}\left({\dot{\omega }}_1\cos \psi -{\dot{\omega }}_2\sin \psi \right)\frac{1}{\sin \theta }\\ -{\dot{\omega }}_1\sin \psi -{\dot{\omega }}_2\cos \psi \\ \frac{\partial \dot{\varphi }}{\partial \psi }\cos \theta \end{array}\right],$$(26) $$\frac{\partial \check{\text{Φ}}}{\partial \dot{\varphi }}=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \dot{\varphi }}+\left[\begin{array}{c}\dot{\theta }\cot \theta \\ 0\\ \dot{\theta }\sin \theta -\frac{\partial \ddot{\varphi }}{\partial \dot{\varphi }}\cos \theta \end{array}\right],$$(27) $$\frac{\partial \ddot{\text{Φ}}}{\partial \dot{\theta }}=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \dot{\theta }}+\left[\begin{array}{c}\frac{1}{\sin \theta }(\dot{\psi }-\dot{\varphi }\cos \varphi )\\ 0\\ \dot{\varphi }\sin \theta -\frac{\partial \ddot{\varphi }}{\partial \dot{\theta }}\cos \theta \end{array}\right],$$(28)

and $$\frac{\partial \ddot{\text{Φ}}}{\partial \dot{\psi }}=\left[\begin{array}{ccc}\frac{\sin \psi }{\sin \theta }& \frac{\cos \psi }{\sin \theta }& 0\\ \cos \psi & -\sin \psi & 0\\ 0& 0& 1\end{array}\right]\frac{\partial \dot{\omega }}{\partial \dot{\psi }}+\left[\begin{array}{c}\frac{1}{\sin \theta }\dot{\theta }\\ -\dot{\varphi }\sin \theta \\ -\frac{\partial \dot{\varphi }}{\partial \dot{\psi }}\cos \theta \end{array}\right],$$(29)

where $\frac{\partial \dot{\omega }}{\partial \varphi },\frac{\partial \dot{\omega }}{\partial \theta }$, and $\frac{\partial \dot{\omega }}{\partial \psi }$ are the row elements of Jacobian $\frac{\partial \dot{\omega }}{\partial \text{Φ}}$, and $\frac{\partial \dot{\omega }}{\partial \dot{\varphi }},\frac{\partial \dot{\omega }}{\partial \dot{\theta }}$, and $\frac{\partial \dot{\omega }}{\partial \dot{\psi }}$ are the row elements of Jacobian $\frac{\partial \dot{\omega }}{\partial \dot{\text{Φ}}}$ respectively.

In order obtain the quantities $\frac{\partial \dot{\omega }}{\partial \text{Φ}}$ and $\frac{\partial \dot{\omega }}{\partial \dot{\text{Φ}}}$, we should again differentiate Eq. (13), and this time with respect to Euler angles and their rates. We get: $$\begin{array}{l}\frac{\partial \dot{\omega }}{\partial \text{Φ}}=\frac{\partial I^{-1}}{\partial \text{Φ}}\left(N-\frac{\text{d}I}{\text{d}t}\omega -\omega \times I\omega \right)\hfill \\ +I^{-1}\left(\frac{\partial N}{\partial \text{Φ}}-\frac{\text{d}I}{\text{d}t}\omega -\omega \times I\omega \right)\hfill \\ +I^{-1}\left(N-\frac{\partial \dot{I}}{\partial \text{Φ}}\omega -\dot{I}\frac{\partial \omega }{\partial \text{Φ}}-\omega \times I\omega \right)\hfill \\ +I^{-1}\left(N-\frac{\text{d}I}{\text{d}t}\omega -\frac{\partial \omega }{\partial \text{Φ}}\times I\omega -\omega \times \frac{\partial I}{\partial \text{Φ}}\omega -\omega \times I\frac{\partial \omega }{\partial \text{Φ}}\right),\hfill \end{array}$$(30)

and $$\begin{array}{l}\frac{\partial \dot{\omega }}{\partial \dot{\text{Φ}}}=I^{-1}\left(-\frac{\partial \dot{I}}{\partial \dot{\text{Φ}}}\omega -\dot{I}\frac{\partial \omega }{\partial \dot{\text{Φ}}}-\omega \times I\omega \right)\hfill \\ +I^{-1}\left(N-\frac{\text{d}I}{\text{d}t}\omega -\frac{\partial \omega }{\partial \dot{\text{Φ}}}\times I\omega -\omega \times I\frac{\partial \omega }{\partial \dot{\text{Φ}}}\right),\hfill \end{array}$$(31)

where $\frac{\partial \omega }{\partial \text{Φ}}$ and $\frac{\partial \omega }{\partial \dot{\text{Φ}}}$ can be formed directly from Eq. (10).

In principle any term in Eqs. (30) and (31) should be modeled as such to calculate the variational equations, practically, we used Maple’s (Char et al. 2013) built-in mathematical algorithms for symbolic computation to calculate partial derivatives of the complex formulas, such as $\frac{\partial \dot{I}}{\partial \text{Φ}}$, and $\frac{\partial I^{-1}}{\partial \text{Φ}}$.

For consistency, with respect to $\frac{\partial I}{\partial \text{Φ}},\frac{\partial \dot{I}}{\partial \text{Φ}},\frac{\partial I^{-1}}{\partial \text{Φ}}$, and $\frac{\partial N}{\partial \text{Φ}}$ expressed in the equations above, we must evaluate the position and velocity components in Phobos’ principle axis system. Here, taking $\frac{\partial I}{\partial \text{Φ}}$ as an example, we have: $$\begin{array}{l}\frac{\partial I}{\partial \text{Φ}}=\frac{\partial I_{\text{tide }}}{\partial \text{Φ}}+\frac{\partial I_{\text{spin }}}{\partial \text{Φ}}\hfill \\ =\frac{\partial I_{\text{tide }}}{\partial \xi }\frac{\partial \xi }{\partial \text{Φ}}+\frac{\partial I_{\text{spin }}}{\partial \text{Φ}}\hfill \\ =\frac{\partial I_{\text{tide }}}{\partial {\xi }_1}\frac{\partial {\xi }_1}{\partial \text{Φ}}+\frac{\partial I_{\text{tide }}}{\partial {\xi }_2}\frac{\partial {\xi }_2}{\partial \text{Φ}}+\frac{\partial I_{\text{tide }}}{\partial {\xi }_3}\frac{\partial {\xi }_3}{\partial \text{Φ}}+\frac{\partial I_{\text{spin }}}{\partial \text{Φ}},\hfill \end{array}$$(32)

where the full descriptions of I_tide and I_spin can refer to Williams et al. (2001). The terms $\frac{\partial \xi }{\partial \text{Φ}}$ represent the change in the principle axis coordinates of a perturbing body (here is Mars) with respect to changes in the Euler angles.

Carrying out the multiplication for the rotation matrices R_ɀ and R_x, we obtain the elements of ℛ, $$\left[\begin{array}{ccc}\cos \psi \cos \varphi -\sin \psi \cos \theta \sin \varphi & \cos \psi \sin \varphi +\sin \psi \cos \theta \cos \varphi & \sin \psi \sin \theta \\ -\sin \psi \cos \varphi -\cos \psi \cos \theta \sin \varphi & -\sin \psi \sin \varphi +\cos \psi \cos \theta \cos \varphi & \cos \psi \sin \theta \\ \sin \theta \sin \varphi & -\sin \theta \sin \varphi & \cos \theta \end{array}\right].$$(33)

Now, r does not depend on the orientation of Phobos’ principle axis coordinates, so: $$\frac{\partial \xi }{\partial \text{Φ}}=\frac{\partial ℛ}{\partial \text{Φ}}r.$$(34)

Thus, we have $$\begin{array}{l}\begin{array}{l}\frac{\partial {\xi }_1}{\partial \varphi }=(-\cos \psi \sin \varphi -\sin \psi \cos \theta \cos \varphi )x\hfill \\ +(\cos \psi \cos \varphi -\sin \psi \cos \theta \sin \varphi )y\hfill \end{array}\\ \begin{array}{l}\frac{\partial {\xi }_1}{\partial \theta }=\sin \psi \sin \theta \sin \varphi x-\sin \psi \sin \theta \cos \varphi y+\sin \varphi \cos \theta z\hfill \\ \frac{\partial {\xi }_1}{\partial \psi }=(-\sin \psi \cos \varphi -\cos \psi \cos \theta \sin \varphi )x\hfill \\ -(\sin \psi \sin \varphi -\cos \psi \cos \theta \cos \varphi )y+\cos \psi \sin \theta z\hfill \\ \frac{\partial {\xi }_2}{\partial \varphi }=(\sin \psi \sin \varphi -\cos \psi \cos \theta \cos \varphi )x\hfill \\ -(\sin \psi \cos \varphi +\cos \psi \cos \theta \sin \varphi )y\hfill \end{array}\\ \begin{array}{l}\frac{\partial {\xi }_2}{\partial \theta }=\sin \psi \sin \theta \sin \varphi x-\sin \psi \sin \theta \cos \varphi y+\sin \varphi \cos \theta z\hfill \\ \frac{\partial {\xi }_2}{\partial \psi }=(-\sin \psi \cos \varphi -\cos \psi \cos \theta \sin \varphi )x\hfill \\ -(\sin \psi \sin \varphi -\cos \psi \cos \theta \cos \varphi )y-\sin \psi \sin \theta z\hfill \\ \frac{\partial {\xi }_3}{\partial \varphi }=\sin \theta \cos \varphi x-\cos \theta \cos \varphi y\hfill \end{array}\\ \begin{array}{l}\frac{\partial {\xi }_3}{\partial \theta }=\cos \theta \sin \varphi x-\cos \theta \sin \varphi y-\sin \theta z\hfill \\ \frac{\partial {\xi }_3}{\partial \psi }=0.\hfill \end{array}\end{array}$$(35)

This concludes the derivation of the variational equations for our Phobos’ rotational model. For the sake of completeness, we have presented most of the important terms in variational equations for rotational equations of motion. By combining the variational equations for orbital equations of motion and geophysical parameters provided in the previous literatures (Lainey et al. 2004b, 2007), we can build the complete variational equations. In this work, only the positions, velocities, Euler angles, and the rate of Phobos are fitted, while their derivatives are integrated simultaneously with the equations of motion and rotation.

4 Comparison of the new dynamical model with the current model

The impetus for the refinements to the Phobos’ rotational and orbital motion described in this paper are the possible availability of high-precision observations from a future mission, such as trajectory data coming from the probe and the image data of Phobos when the probe is on very close orbits. In addition, a lander and tracking measurement carried out on Phobos would prove especially valuable (Kawakatsu et al. 2017; Usui et al. 2018). In order to verify the feasibility and effectiveness of the newly constructed model, we need to distinguish the level at which current models fail as a result of not considering the fully 3D rotation of Phobos.

4.1 Reference ephemerides

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. 2021), 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 0 = (2e + 𝓐) sin M, where e, 𝓐 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: $$F_m=\frac{3}{2}{\mu }_p\left(\frac{a_p^2}{r^4}\right)\left[\left(-C_{20}+6C_{22}\cos 2\theta \right)\widehat{r}+4C_{22}\sin 2\theta \widehat{t}\right],$$(36)

where µ_p is the GM of Phobos, a_p is the equatorial radius of Phobos, C₂₀ and C₂₂ are the second degree harmonic of Phobos, $\widehat{r}$ is the unit vector directed from Mars toward satellite, and $\widehat{t}$ is the unit vector in the satellite’s orbit plane normal to $\widehat{r}$ and in the direction of its orbital motion. Easily, the reactive force acting on the satellite is: $$F_s=-\left(\frac{{\mu }_m}{{\mu }_p}\right)F_m,$$(37)

with µ_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. 2021) as the “observations” and fit our model to them. The adjustment was carried out in Cartesian planeto-centric 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, Si, the C₂₀ and C₂₂ 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.

4.2 Adjustment to the simpler model

To explore the difference between our full 3D dynamical model and the previous simple model (Jacobson 2010; Jacobson & Lainey 2014; Lainey et al. 2021). We adjusted the full rotational model to the reference ephemeris integrated in Sect. 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.

Fig. 1 Difference in distance after fitting the numerical model to the NOE-4-2020 ephemerides for Phobos.

Fig. 2 Residuals after an adjustment over 10 yr of our full model to the simulated simpler model.

Fig. 3 Difference in distance after 10 yr of fitting the full model to the simulated simpler model.

Fig. 4 Temporal evolution of Phobos’ Euler angles about 10 yr.

Fig. 5 Angles between the direction from Phobos to Mars and Phobos’s axis of minimum principal moment of inertia.

Fig. 6 Angles between the pole of Phobos and normal of Phobos’s orbital plane.

Fig. 7 Temporal evolution of Phobos’ pole position on inertial space about 10 years.

4.3 Full model with higher gravity

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₃₃ to be the cause of large orbit difference, that is, if we neglect the S₃₃ and the adjustment will be very close to the case where only the second-degree gravity field is considered. Figure 9 presents the difference in distance after fitting the full model with the third-degree gravity field, but not S₃₃ 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₃₃, and the tidal orbital expansion may also bring about this signal (Lainey et al. 2020). 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, ⋯ , 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.

Fig. 8 Difference in distance after 10 years of fitting the full model with gravity filed truncated at the third degree to the simulated simpler model.

Fig. 9 Difference in distance after 10 yr of fitting the full model with gravity filed truncated at the third degree but without S33 to the simulated simpler model.

Fig. 10 Difference in distance after 10 yr of fitting the full model with gravity filed truncated at the fourth and third degrees to the simulated simpler model. Note: the unit of the Y-axis is meters.

4.4 Effect of Phobos’ tidal Love number, k₂

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₂, was referred to by Le Maistre et al. (2013) to construct the numerical model of motion. The value of k₂ 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₂ on Phobos’ orbital motion, we chose various k₂ to simulate Phobos’ motion.

By including a range of plausible value for k₂ ∊ [1.83 × 10⁻⁷, 5.3 × 10⁻⁴] (Le Maistre et al. 2013), our new model uses a second-degree gravity field of Phobos and a different k₂ value fitted to the simulated simple model that without considering k₂ 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 m in ten years’ integration when k₂ ≤ 1.0 × 10⁻⁴. This result suggests that it is difficult to effectively determine the k₂ value of Phobos from the current data.

Fig. 11 Difference between the integration results of the full model using the higher order gravity field and the fourth-degree model, with the initial values all using the fitted values of the fourth-degree model. Note: the unit of the Y-axis is meters.

Fig. 12 Post-fitted orbital position differences over 10 years between using different k2 with k2 = 1.0 × 10−4. Note: the unit of the Y-axis is meters.

5 Discussion and conclusion

A high-precision numerical ephemeris can typically provide the accurate positions, velocities, and orientations of celestial bodies over time. Therefore, it can allow for the detailed study of the evolution and internal structure of these celestial bodies. In this paper, we establish a numerical dynamical model of Phobos’ motion with full consideration of its rotation. The full “3D” dynamical model was constructed in planetocentric Cartesian coordinates based on integrating Newton’s equations for orbital motion and Euler’s rotational equations for Phobos’ rotation simultaneously. We explicitly gave the variational equations of the rotational motion of Phobos when adjusting the numerical solution to the observations. In order to exclude the influence of other parameters in the model, we first simulated the simple model currently in use and then fitted it to the current French ephemerides NOE-4-2020, using the least-squares approach to determine the initial conditions.

We fit full rotation model to the simulated simple model, with the gravity field truncated at only C₂₀ and C₂₂. The results of the simple and full models are close, but the full model includes more latitudinal motion in the rotational motion. However, when we fit the full model including the third-order gravity field coefficients to the simulated simple model, the results show a large discrepancy, this may be due to the introduction of new librations resulting from the use of gravity field coefficients containing higher orders, a phenomenon that needs to be re-examined in the future based on more plausible gravity field coefficients.

Given the internal structure and elastic property of Phobos have not been clearly determined so far, we introduced the tidal Love number, k₂, of Phobos into our model to check its influence on Phobos’ motion. The simulations indicate that this perturbation is non-negligible when the k₂ value is larger than 1.0 × 10⁻⁴. Otherwise, the satellite can be treated as a rigid body when the k₂ is no larger than 1.0 × 10⁻⁴, which will make the integration and fit process more efficient and easier.

This revised numerical model of Phobos’ motion as an extension of the previous work of Lainey et al. (2007); Jacobson (2010); Jacobson & Lainey (2014) and it provides potential opportunity for future use of high-precision observation from future missions to further the study of Phobos and Mars. In the future, it is not only the position and orientation that will feasibly be determined by using this new model, but also Phobos’ gravity coefficients and k₂, for instance, when combined with observation data from future mission. Finally, it needs to be pointed out is that our modeling process is based on universal principles and can be applied to a full Mars system (including Deimos) or other planetary systems to confront the continuous improvement of the accuracy of observation data.

Acknowledgements

This research is supported by the National Key Research and Development Program of China (2021YFA0715101, 2022YFF0503202), the National Natural Science Foundation of China (12033009, 12103087, 42241116), the International Partnership Program of Chinese Academy of Sciences (020GJHZ2022034FN), the Yunnan Province Foundation (202201AU070225, 202301AT070328), the Young Talent Project of Yunnan Revitalization Talent Support Program, grant from Key Laboratory of Lunar and Deep Space Exploration, Chinese Academy of Sciences (LDSE202004), and open research fund of state key laboratory of information engineering in surveying, mapping and remote sensing, Wuhan University (21P02). Y. Yang thanks Dr. Valéry Lainey for the careful reviewing and nice suggestions that have been incorporated into and improved the paper and kindly provided the subroutine that compute Phobos’ ephemeris from NOE ephemerides, thank the discussion with Dr. Jinling Li (about computing the variational equations of rotational motion). Phobos’ ephemerides files can be downloaded from ftp://ftp.imcce.fr/pub/ephem/NOE/(NOE-4-2020).

References

All Tables

All Figures

	Fig. 1 Difference in distance after fitting the numerical model to the NOE-4-2020 ephemerides for Phobos.
In the text

	Fig. 2 Residuals after an adjustment over 10 yr of our full model to the simulated simpler model.
In the text

	Fig. 3 Difference in distance after 10 yr of fitting the full model to the simulated simpler model.
In the text

	Fig. 4 Temporal evolution of Phobos’ Euler angles about 10 yr.
In the text

	Fig. 5 Angles between the direction from Phobos to Mars and Phobos’s axis of minimum principal moment of inertia.
In the text

	Fig. 6 Angles between the pole of Phobos and normal of Phobos’s orbital plane.
In the text

	Fig. 7 Temporal evolution of Phobos’ pole position on inertial space about 10 years.
In the text

	Fig. 8 Difference in distance after 10 years of fitting the full model with gravity filed truncated at the third degree to the simulated simpler model.
In the text

	Fig. 9 Difference in distance after 10 yr of fitting the full model with gravity filed truncated at the third degree but without S33 to the simulated simpler model.
In the text

	Fig. 10 Difference in distance after 10 yr of fitting the full model with gravity filed truncated at the fourth and third degrees to the simulated simpler model. Note: the unit of the Y-axis is meters.
In the text

	Fig. 11 Difference between the integration results of the full model using the higher order gravity field and the fourth-degree model, with the initial values all using the fitted values of the fourth-degree model. Note: the unit of the Y-axis is meters.
In the text

	Fig. 12 Post-fitted orbital position differences over 10 years between using different k2 with k2 = 1.0 × 10−4. Note: the unit of the Y-axis is meters.
In the text