1. Introduction
Human understanding of the Martian moons began with the discovery by Asaph Hall, a scientist at the United States Naval Observatory, in 1877. He identified the two moons of Mars and named them Phobos and Deimos, respectively. Researchers have conducted comprehensive studies on Phobos [
1,
2]. Compared to Phobos, there is relatively limited publicly available research on Deimos. Therefore, this paper focuses on Deimos as the primary research target. The physical parameters of Deimos are shown in
Table 1 [
3]. The study of the Martian moons not only contributes to a deeper understanding of the formation and evolution of terrestrial planet systems but also provides crucial clues in researching the formation and evolution of the solar system. Up to now, humans have launched multiple probes to Mars and conducted extensive exploration activities on its two moons. Approximately 90 years after the initial discovery of the two Martian moons, the first artificial satellite dedicated to Mars, Mariner 4, conducted a flyby of Mars in 1965 [
4]. Subsequently, human exploration of Mars has not ceased, and spacecraft have been launched to Mars on multiple occasions. Mariners 6 and 7 also carried out flyby missions of Mars. In 1971, the first Mars orbiter, Mariner 9, was launched. It obtained 214 images of Phobos and Deimos during close flybys, marking the earliest spacecraft observations of the Martian moons. In the late 1970s, the Soviet Union successively launched two orbiters, Viking 1 and Viking 2, to Mars. These spacecraft were equipped with landers. Among them, Mars Express, launched by the European Space Agency (ESA), acquired the most data on the Martian moons. Its impact on understanding the orbits of the Martian moons and the Martian ephemerides has been profound. Notably, data from the Pathfinder Mars lander played a crucial role in the study of Mars’ rotation and gravity field. The gravity field model MRO120F from the Pathfinder mission is utilized in this research [
5]. In recent years, both China and Japan have initiated their own Mars exploration missions [
6,
7,
8]. It is anticipated that an abundance of high-precision data will be instrumental in advancing the research on the dynamical models of Martian moons.
In various exploration activities of celestial bodies beyond Earth, investigating their orbital trajectories is a crucial scientific inquiry, holding significant academic importance for a profound understanding of the internal structure and physical parameters of these celestial entities. The data collected from ground-based observations to spacecraft missions provide the foundation for our research on the dynamical models of Martian moons. Multiple international institutions have conducted in-depth research in constructing dynamical models for Martian moons. With the improvement in the volume and precision of human exploration data for Martian moons, the dynamical models for these moons have undergone a developmental process from analytical models and semi-analytical models to numerical models.
In the early stages, dynamical models for Martian moons were typically established using analytical methods. The earliest analytical model was developed by Struve in 1911. He fitted observational data from the United States Naval Observatory collected between 1877 and 1909, calculating numerical values for the Martian moon’s mean longitude, orbital elements, mean motion, and precession rate. It was not until Shor that an analytical theory was introduced that incorporated the influence of tidal acceleration into the model [
9]. This modification involved adding an acceleration term to the mean orbital longitude, and fitting was conducted using observational data spanning from 1877 to 1973.
With the continuous improvement in observational precision, analytical methods gradually revealed their limitations, leading to the brief emergence of semi-analytical approaches in dynamical modeling. The theory proposed by Chapront-Touzé [
10,
11,
12] represents a high-precision semi-analytical model. It incorporates perturbations from the Sun and planets, the gravitational field of Mars, and mutual perturbations between Phobos and Deimos. This model is a more comprehensive semi-analytical approach to constructing dynamical models for Martian moons. Through fitting observational data from both ground-based and spacecraft observations spanning from 1877 to 1988, an ephemeris for the Martian moons was established.
With the development of deep space exploration, especially the increase in observational data and the continuous improvement in computing power, numerical methods, characterized by high computational accuracy and simplicity of formulas, have gradually replaced analytical methods as the primary approach in planetary orbit calculations. In 1977, Tolson [
13] used numerical methods and Viking spacecraft data to generate the first ephemeris for Phobos. In 2007, Lainey et al. [
14] developed the first complete numerical integration ephemeris for the Martian moon. The dynamic model included a 10th-order non-spherical Martian gravity field [
15], perturbations from the Sun, Jupiter, Saturn, Earth, and the Moon, mutual perturbations between the moons, tidal effects caused by the two moons on Mars, and the quadrupole moment of Phobos. Subsequently, Jacobson [
16,
17] further incorporated the librational effects of the moons in the dynamical model, considering perturbations caused by the oscillation of the axis pointing toward Mars in the moon’s body-fixed coordinate system within the orbital plane. Additionally, a new Mars rotation model was adopted, and the model was fitted with the Mars Express dataset, generating the Martian moon ephemeris MARS097.
In Jacobson’s dynamical model, while the libration of the minimum moment of principal axis of Deimos in the orbital plane is considered, it does not incorporate the full rotation of the moon as is the case in the lunar dynamical model. Therefore, the aim of this paper is to establish a rotational model for Deimos in the inertial frame based on the Euler–Liouville rotational equations, following the approach used in the lunar dynamical model. The purpose is to further develop the simple dynamical model for Martian moons and establish a full dynamical model for Deimos in the inertial frame. Simultaneously, based on the rotation model for Deimos, the study incorporates the influence of Deimos’ own gravitational field when it orbits Mars, completing the full dynamical model of the rotational and orbital coupling for Deimos.
The first section of this paper will introduce the establishment process and distinctions between the simple model and the dynamical model considering the coupling of rotation and orbit. In the second section, the least squares method will be employed to introduce the data-fitting generation process and the adjustment model for precise orbit determination using the established dynamical model. Finally, a comparison will be made between the changes in coordinate axes during the orbital integration in inertial space for the two models.
2. The Establishment of the Simple Model and Full Model
2.1. The Rotation of the Mars
When describing the motion of Martian moons, the coordinate system chosen is the International Celestial Reference Frame (ICRF), while the gravitational field is represented in the Martian body-fixed coordinate system. Therefore, calculating the acceleration induced by the Martian gravity field on the moons requires an appropriate rotation model. Konopliv et al. conducted an in-depth study of the Mars rotation based on data from the Pathfinder lander, establishing the most accurate current rotation model for Mars [
5]. According to this model, the method for converting positions from the Mars-fixed frame to the ICRF is as follows:
where
and
represent the coordinates of Deimos in the inertial frame and the body-fixed frame, respectively.
N is the angle in the plane of the Earth-mean-equator of J2000 (EME2000) from the vernal equinox (intersection of the mean ecliptic and EME2000 planes) to the intersection of the Mars-mean-orbit of J2000 and EME2000 planes;
J is the inclination of the Mars-mean-orbit plane relative to the EME2000 plane;
is the angle in the Mars-mean-orbit plane from the intersection of the EME2000 and Mars-mean-orbit planes to the intersection of the Mars-mean-orbit and Mars-true-equator of date planes;
I is the inclination of the Mars true equator relative to the Mars-mean-orbit plane;
is the angle in the plane of the Mars true equator from the intersection of the Mars-mean-orbit plane and Mars true equator to the Mars prime meridian; and
,
represents the polar motion.
2.2. Mars Gravity Field
Due to the irregular shape and uneven mass distribution of Mars, when calculating the gravity on Deimos, it cannot be treated as an ideal homogeneous sphere. Non-spherical perturbations should be considered. The gravitational potential function is expanded into spherical harmonics as follows:
In this expression, is the gravitational constant of Mars, is the radius of Mars, represents the normalized associated Legendre function, and and are the normalized spherical harmonic coefficients of the Mars gravity field model. , , and r denote the longitude, latitude, and distance to the center of Mars in the fixed coordinate system of Deimos on Mars. When using the gravity field formula, it is necessary to follow the coordinate transformation method described in the previous section. First, the coordinates of Deimos in the inertial frame should be converted to the Mars-fixed frame. Subsequently, the integration calculations are performed by transforming from the Mars-fixed frame back to the inertial frame.
When computing gravitational field acceleration, a recursive method is commonly employed to avoid complex differentiation processes. For the Deimos orbit calculation, taking the example of a 10th-order Martian gravity field, the Cunningham [
18] algorithm is used to calculate the gravitational field and its partial derivatives.
When calculating the Martian gravity potential, Cunningham provided a single recursion formula that can be used to calculate the gravity potential and the corresponding acceleration in Cartesian coordinates. The formula is defined as follows:
where the gravitational field can be represented as (non-normalized coefficients are used here):
Then, the acceleration of the gravitational field can be represented as:
2.3. General Relativity Perturbation Model
During the motion of Deimos, it is mainly influenced by the General Relativity effect generated by the mass of Mars. Since the eccentricity of Deimos’ orbit is very small, it can be considered a nearly circular orbit. In this case, the acceleration induced by the General Relativity effect of Mars on Deimos is given by:
where
is the gravitational constant of Mars,
v is the velocity of Deimos in the Mars inertial frame,
is the position vector from the Mars to Deimos, and
c is the speed of light.
2.4. Solid Tide Perturbation
Due to the influence of the gravitational forces from other celestial bodies in the solar system, Mars experiences changes in its gravitational field, causing solid tidal effects. For a celestial body with a mass
and a radius
, the solid tidal effect (according to elastic model pointed in [
19] and also [
20]) on Deimos’ gravitational on Deimos’ gravitational potential function is given by [
21]:
where
is the position vector from the perturbing body to Deimos,
is the position vector from Mars to Deimos,
is the position vector from Mars to the perturbing body, and the Love number of Mars has a value of 0.169 [
20]. The solid tide model used in this study does not incorporate the tidal energy dissipation effects of Mars. We referred to previous studies on tidal energy dissipation [
22,
23,
24,
25,
26], and integrating the established dynamical model over ten years in numerical experiments reveals that tidal energy dissipation effects only cause differences of less than a few meters. Therefore, within the time scale of this study (10 years), we have not considered the influence of tidal dissipation.
2.5. Simple Model
The current Deimos dynamical model is referred to as a simplified model. In the simple Deimos dynamical model, the librations of Deimos are considered. Jacobson [
16] introduced the specific expression for this effect:
where
represents the gravitational constant of Mars;
R is the radius of Deimos;
r represents the distance from Mars to Deimos;
is the second-degree sectoral harmonic coefficient of Deimos;
is the second-degree tesseral harmonic coefficient of Deimos;
f and
M are the true anomaly and mean anomaly of Deimos’ orbit, respectively;
e is the eccentricity of the orbit;
A,
B, and
C are the moments of inertia with respect to the principal axes of Deimos;
is the unit vector from Mars to Deimos; and
is the unit vector in the Deimos orbital plane that is perpendicular to the velocity vector of Deimos. This model is the dynamical model used in the current ephemeris.
2.6. Full Model
The main difference between the full model and the simple model lies in the consideration of the complete rotation of Deimos in inertial space in the full model. There have been some research efforts on the rotation of Phobos [
27,
28]. In this section, we will establish the rotation model of Deimos based on the theory of rigid body rotation. Based on the established rotation model, it is possible to calculate the gravitational force of Deimos on Mars and the acceleration of Deimos caused by the reaction force it experiences due to the gravitational interaction with Mars. Thus, a dynamical model for the coupled rotation and orbit of Deimos has been established. The gravitational field coefficients for Deimos are listed in
Table 2 [
3]. This study further advances the current dynamical model that only considers the simple librations of Deimos. It comprehensively calculates the rotation of Deimos in inertial space due to torques from celestial bodies such as Mars and the Sun.
The differential equations governing the rotational motion of a rigid body are classically represented by the Euler–Liouville equations, and their form in an inertial frame is as follows:
where
N is the torque acting on Deimos,
represents the rotation of Deimos, and
I is inertia tensor matrix.
When describing the rotation of Deimos, we take the center of mass of Deimos as the origin, the three principal axes of Deimos as the coordinate axes, and establish the Deimos Principal Axis (PA) reference frame rotating with Deimos. Simultaneously, we can obtain the celestial reference frame of the Deimos center of mass with coordinates aligned with ICRS. Assuming Deimos as an ideal rigid body, the transformation between the Deimos body-fixed coordinate system and the inertial frame is determined by three Euler angles (
,
,
) describing the rotation of the rigid body. This transformation process is expressed as:
represents the rotation of the rigid-body moon under the action of torque in the inertial coordinate system, and the relationship between
and the three Euler angles and the rates of change of the three Euler angles can be obtained:
Solving for the first derivative of
, we can obtain the accelerations of the three Euler angles:
However, when employing Equation (
16), one may encounter singular situations. In such cases, using Wisdom angles instead of Euler angles might be more beneficial [
29]. According to the Euler–Liouville equation, the value of
depends on the torque acting on the Deimos and the moment of inertia. Taking the first derivative of Equation (
13) yields the expression:
I is the inertia tensor matrix, which can be expressed as:
The torque acting on Deimos during its rotation can be divided into three components: (1) the torque produced by Mars as a point mass in the Deimos gravitational field; (2) the torque exerted by the Sun, treating it as a point mass due to the presence of the Deimos gravitational field; and (3) the torque generated in the Deimos due to the irregular shape of Mars, under the influence of the Deimos gravitational field [
30]. The formulas can be expressed as:
where
and
represent the torques generated by the point masses Mars and the Sun, and
represents the torque produced by the shape effect of Mars [
31].
is generated by the second-order term of the Martian gravitational field and is mainly related to the value of Mars. In numerical experiments, the torque generated by Phobos is a magnitude smaller than that of Mars and the Sun. Therefore, the model does not consider the impact of the torque generated by Phobos. The force exerted by Mars in the moon’s gravitational field is expressed as:
The torque produced by this force in the Deimos is expressed as:
represents the mass of Mars,
represents the coordinates of Deimos in the Mars-centered inertial frame, and
represents the gravity field of Deimos. And
represents the torque induced by the oblateness of Mars and the shape effect of Deimos. The expression is given by [
30]:
is the mass of Deimos, represents the value of the Mars flattening, is the distance between Deimos and Mars, is the unit vector pointing from Deimos to Mars, is the direction of Mars’s polar axis, and all vectors are described in the Deimos body-fixed coordinate system.
According to the established model, the initial Euler angles of Deimos at the J2000 epoch are provided by Brent et al. [
32]. The integration of the established rotational differential equations using the ABM integration algorithm yields the results for the three Euler angles as shown in
Figure 1 and
Figure 2.
2.7. Result of the Model
At this point, we have considered factors influencing the motion of Deimos, utilizing a second-order gravity field for Deimos. We have established a dynamical model for the coupled rotation and orbit of Deimos. The initial value of Deimos’ position at the J2000 epoch is provided by MARS097. The integration orbit for both the simple dynamical model and the full model with rotation–orbit coupling established in this study are compared in
Figure 3.
The experimental results indicate that, considering the J2000 epoch as the integration starting point with a time span of 10 years, the difference in the computed results between the two dynamical models is around 10 km. The difference increases linearly over time, primarily due to the inclusion of rotation and the gravitational field of Deimos itself in the dynamical model. This discrepancy suggests that considering the complete rotation effect of Deimos in the dynamical model is essential for improving its accuracy.
When considering the third-order gravity field of Deimos, the integrated orbit of Deimos shows differences compared to the integrated orbit of the dynamical model considering only the second-order gravity field of Deimos as shown in
Figure 4.
The model results indicate that, after integrating for 10 years, the difference between the second- and third-order gravity fields of Deimos is only on the order of meters, which is an extremely small value. Therefore, when constructing a comprehensive model for Deimos, considering gravity field orders up to the second order is sufficient.
4. Conclusions
This study presents a novel model describing the motion of Deimos in an inertial frame. Building upon existing models, it thoroughly calculates the rotational effects induced by the irregular shape of Deimos itself and the torques exerted by large celestial bodies (Mars, Sun) in the solar system in inertial space. Furthermore, it considers the coupling effect resulting from rotation and orbit, establishing a full dynamical model for the rotation of Deimos coupled with its orbit. The model is integrated using the Adams–Bashforth–Moulton (ABM) integration algorithm to concurrently solve the 12 elements governing rotation and orbital dynamics in the differential equations.
To validate the model, precise orbit determination techniques for artificial satellites are employed to fit data to the full model. The results demonstrate the stability and reliability of the established model. Additionally, the study derives an analytical expression for the solution of the coupling effect between Deimos’ rotation and orbit. The variational equation, starting from the J2000 epoch, is integrated over ten years. The final outcomes reveal differences of less than 50 m when compared to the current Deimos ephemeris MARS097, showcasing its practical applicability in engineering.
As various exploration missions targeting the Martian moons continue in the future, it is anticipated that a considerable amount of high-precision data on these moons will be generated. Leveraging these data in conjunction with the full dynamical model established in this study, there is potential to create a new generation of more accurate ephemerides for the Martian moons. Additionally, by utilizing the full model established in this paper along with future observational data, there is the prospect of quantitatively studying orbital resonances among celestial bodies in the solar system [
34,
35] and the Lense–Thirring effect [
36,
37,
38] in orbital mechanics. This, in turn, holds the potential to further enhance the accuracy of the Deimos’ dynamical model.
The model established in this paper lays the foundation for further constructing an accurate ephemeris for Deimos. It demonstrates versatility in solving the dynamical equations for celestial bodies similar to Deimos and can serve as a reference for developing numerical ephemerides for major celestial bodies in the solar system.