© The Authors 2023

1 Introduction

Apart from the Earth and the Moon, the rotation of Solar System bodies is inferred from data that are analyzed using the barycentric celestial reference system (BCRS). Most of the time, the inferred rotation model is then used to estimate physical properties of these bodies related to their atmosphere and surface dynamics or to their interior structure and composition. To prevent errors in the physical interpretation of the rotation models, it is necessary to properly correct them for relativistic contributions.

Thanks to the many missions that have visited and studied Mars over recent decades, the rotation of the red planet has been thoroughly investigated. Many solutions have been proposed in the literature based on different sets of data (mostly radiometric data from orbiters and landed spacecraft, e.g., Konopliv et al. 2020; Kahan et al. 2021). Along with these rotation models, interior models and global circulation models (GCM) have also been produced and improved over time, especially following the InSight mission. This is why we focus the present paper on Mars, although the levels of relativistic contributions in the rotation of other planets are also provided at the end of this paper.

The orientation of Mars with respect to its orbit can be described with a set of Euler angles (see Fig. 1), which are used to build the rotation models in the BCRS. These models are commonly corrected for relativistic effects in the longitude angle (ψ) and in the spin angle (ϕ), while no relativistic correction in the obliquity (ε) of Mars is usually applied in the literature. The relativistic correction in ψ does not often explicitly appear in the angle definition¹, but instead it is distributed among the different terms used to express that angle (e.g., Konopliv et al. 2006, Eq. (14); Folkner et al. 1997a, Eq. (5); or Baland et al. 2020, Eq. (14)): $$\psi \left(t\right)={\psi }_0+{\dot{\psi }}_{0}t+{\psi }_{\text{nut}}\left(t\right),$$(1)

where ψ₀ is the constant value at epoch J2000, ${\dot{\psi }}_0$ is the constant Mars precession rate, and ψ_nut(t) is a periodic series of nutations. Baland et al. (2020) found a geodetic precession rate of 6.7 mas yr⁻¹ hidden in ${\dot{\psi }}_0$ (see their Eq. (69)), which is greater than the uncertainty on the determination of the precession rate (~2mas yr⁻¹, Le Maistre et al., in prep.), as well as periodic geodetic nutations in ψ_nut, including an annual term with 0.6 mas in amplitude (see also Eroshkin & Pashkevich 2007).

As in the case of the Earth, Mars experiences periodic variations in its rotation due to atmosphere and surface dynamics that result in length-of-day (LOD) variations. By convention, the rotation angle of Mars, ϕ, is measured from the ascending node of the Mars true equator of date over the Mars mean orbit of epoch (usually chosen as the 1980 orbit, following Folkner et al. 1997b) to the intersection of Mars's Prime Meridian on Mars's true equator of date (see Fig. 1). It is generally decomposed as (e.g., Konopliv et al. 2006, Eq. (16)): $$\begin{array}{c}\varphi \left(t\right)={\varphi }_0+{\dot{\varphi }}_{0}t-{\psi }_{nut}\left(t\right)\cos {\epsilon }_0\\ +{\displaystyle \sum _{j=1}^4\left({\varphi }_{cj}\cos \text{\hspace{0.17em}}j{l}^{\prime }\left(t\right)+{\varphi }_{sj}\sin j{l}^{\prime }\left(t\right)\right)}+{\left[\varphi \right]}_{\text{GR}}\left(t\right),\end{array}$$(2)

where ϕ₀ is the constant value at epoch J2000, ${\dot{\varphi }}_0$ is the constant spin rate of Mars, and ϕ_cj and ϕ_sj are the amplitudes of the periodic variations induced by the seasonal atmosphere and surface dynamics, with l′(t) as the mean anomaly for Mars. Additionally, ψ_nut(t) is the periodic nutation in longitude of Eq. (1) and ε₀ is the J2000 epoch Mars obliquity, t is barycentric dynamical time (TDB), as detailed in Sect. 2.2. The term −ψ_nut(t) cos ε₀ in Eq. (2) is a correction for the seasonal periodic variations that are measured along Mars's mean equator of date – rather than along the true equator of date, which is nutating with respect to the mean equator of date. This nutation term in the rotation angle therefore includes the periodic relativistic correction in longitude. The last term [ϕ]_GR(t) in Eq. (2) is another relativistic correction arising when changing the time scale from the proper time of Mars to TDB. This term takes into account: (i) the time dilation effect due to the barycentric orbital velocity of Mars's center-of-mass and (ii) the Einstein gravitational redshift effect due to the change of Mars's altitude inside the gravitational potential of the Sun caused by the non-vanishing eccentricity. This relativistic correction is of the same order of magnitude as the seasonal periodic terms, and must be estimated accurately.

Yoder & Standisti (1997) provided an estimation for [ϕ]_GR(t) (for details, see their Eq. (21), with amplitudes in mas): $$\begin{array}{cc}{\left[\varphi \right]}_{\text{GR}}\left(t\right)\approx {\displaystyle \sum _{j=1}^3{\varphi }_{rj}\sin j{l}^{\prime }\left(t\right)}=-175.80\text{\hspace{0.17em}}\sin l^{\prime }\left(t\right)& \\ -8.20\sin 2l^{\prime }\left(t\right)-0.60\sin 3l^{\prime }\left(t\right),& \end{array}$$(3)

which is still in use (e.g., in Konopliv et al. 2020). This estimation is expressed as a sum of trigonometric terms, the arguments being the harmonics of the Mars mean anomaly, l′(t), as well as for the seasonal terms. As we show in Sect. 3, this estimation is affected by an error of about 9 mas on the periodic terms, which is greater than the current formal uncertainty on the determination of rotation periodic variations (~1 mas, Le Maistre et al., in prep.) and lacks a linear term that would affect the rotation rate.

To ensure a correct interpretation of rotation variations measurements, we aim to estimate the relativistic contribution to the Euler angles at the 0.1 mas level. We go on to update the estimation of the main terms and also investigate the existence of terms neglected thus far with amplitudes higher than the 0.1 mas threshold.

The paper is organized as follows. In Sect. 2, we present the theory for the geodetic and Lense-Thirring effects, as well as for the effect of the time transformation. The results for the geodetic and Lense-Thirring effects are also presented in Sect. 2, whereas the results for the effect of the time transformation are presented in a dedicated section (Sect. 3), as it requires more investigation. In that section, we also present solutions obtained initially based on the assumption that the planets evolve on Keplerian orbits and then based on realistic orbits as described by ephemerides. In Sect. 4, we discuss the signature of the relativistic effects on the Doppler signal of a Martian radioscience instrument. In Sect. 5, we briefly introduce a relevant application of the model to other planets in the Solar System. Our discussion and conclusion are given in Sect. 6.

Fig. 1 Euler angles between the rotating body frame of Mars (axes XYZ) and the inertial frame associated with the mean orbit of Mars of 1980 (axes xyz, the x-axis is in the direction of the ascending node of Mars's orbit over the Earth ecliptic of epoch J2000). The X-axis of the BF is chosen as the prime meridian defined in the IAU convention (Archinal et al. 2018). The spin axis longitude ψ is measured from the x-axis to the autumn equinox, ϕ is measured from the equinox to the X-axis, and the obliquity, ε, is the angle from the z-axis to the Z-axis – or the inclination of the BF equator over the IF xy plane.

2 Theory

Two different types of relativistic contributions exist which arise in the BCRS rotation model of Solar System bodies: (i) contributions that impact directly the rotation of the body and (ii) contributions arising from reference frame transformation. The first type of contributions concerns the ones that impacts the spin equation of motion such as the geodetic precession and nutations (Fukushima 1991; Eroshkin & Pashkevich 2007; Baland et al. 2020) and the Lense-Thirring and Thomas precessions (see Sect. 2.1). The second relativistic contributions come from the reference frame transformation between a local inertial frame that would be used to describe the local physics of the body and the BCRS used to analyze the data. The theory of reference frame transformation to first post-Newtonian order has been derived by Brumberg & Kopejkin (1989); Kopejkin (1988); Damour et al. (1991); Klioner & Voinov (1993) – and it has also been adopted in the IAU 2000 conventions, as in Soffel et al. (2003). Of prime importance for our purpose is the time transformation between a local reference frame and BCRS. In Sect. 2.2, we will present in details various contributions that arise in the time transformation and their impact in terms of the Mars rotation model.

2.1 Geodetic, Lense-Thirring, and Thomas precession effects

Within general relativity framework, the evolution of a spinning body is given, at the first post-Newtonian approximation, by a simple precession relation (see e.g., Barker & O'Connell 1970; Soffel et al. 2003; Poisson & Will 2014): $$\frac{\text{d}S}{\text{d}t}=\Omega \wedge S+\mathcal{O}\left(c^{-4}\right).$$(4)

In this equation, S denotes the spin angular momentum of Mars and Ω is the relativistic total precessional angular velocity, which is decomposed into three parts: Ω = Ω_so + Ω_ss + Ω_ΤΡ. The term Ω_so is called the spin-orbit precessional angular velocity, is the spin-spin precessional angular velocity, and Ω_ΤΡ is called the angular velocity of the Thomas precession. The spin-orbit and spin-spin components are also called the "geodetic precession" and the "Lense-Thirring precession," respectively.

When applying Eq. (4) to the description of the spin variations of Mars, we can keep the contribution of the Sun only in the relations of the spin-orbit and spin-spin precessional angular velocities, while the contribution of Phobos alone can be retained in the relation of the Thomas precession angular velocity. Thus, the expressions of the angular velocities are as follows: $${\Omega }_{\text{so}}=\frac{G{M}_{\odot }}{2c^2{\parallel X-X_{\odot }\parallel }^3}\left(X-X_{\odot }\right)\wedge \left(3v-4v_{\odot }\right),$$(5a) $${\Omega }_{\text{ss}}=\frac{G{S}_{\odot }}{c^2{\parallel X-X_{\odot }\parallel }^3}\left(\frac{3\left(\left(X-X_{\odot }\right)\cdot {\widehat{e}}_{\odot }\right)\left(X-X_{\odot }\right)}{\parallel X-X_{\odot }{\parallel }^3}-{\widehat{e}}_{\odot }\right)$$(5b) $${\Omega }_{\text{TP}}=-\frac{1}{2c^2}\text{v\hspace{0.17em}}\wedge Q,$$(5c)

where x and x_⊙ are the barycentric positions of Mars and of the Sun, respectively, and v and v_⊙ represent their barycentric velocities. The vector Q is the non-geodesic acceleration whose expression can be derived from Eq. (6.30a) of Damour et al. (1991); hereafter, we only consider the dominant Newtonian contribution. In addition, M_⊙ is the mass of the Sun, while S_⊙ and ê_⊙ denote the magnitude and direction of the Sun's spin angular momentum, respectively, c is the speed of light, and G is the universal gravitational constant.

Hereafter, we write ê_⊙ = cos α_⊙ cos δ_⊙ ê_x + sin α_⊙ cos δ_⊙ ê_y + sin δ_⊙ ê_z with α_⊙ and δ_⊙ as the right ascension and declination of the direction of the Sun's spin axis, respectively. Both are assumed to be fixed in the inertial frame associated with the mean orbit of Mars, namely (ê_x, ê_y, ê_z). We also consider that M/M_⊙ ≪ 1, so that Mars follows an heliocentric orbit, namely ‖x_⊙‖/‖x‖ ≪ 1, ‖v_⊙‖/‖v‖ ≪ 1. We denote by ê_z the direction of the orbital angular momentum of Mars (assumed to be constant) and by ê_x the direction of the ascending node of Mars's orbit in the ecliptic; ê_y completes the triad such that (ê_x, ê_y, ê_z) is a direct orthogonal basis. After decomposing Eqs. (5a) and (5b) into this basis (the case of the Thomas precession is treated separately at the end of the section), we obtain the following relationships: $${\Omega }_{\text{so}}\text{}={\text{Ω}}_{\text{so}}{\left(1+e\text{\hspace{0.17em}}\cos f\right)}^3{\widehat{e}}_{\text{z}}\text{}={\text{Ω}}_{s\text{o}}^z{\widehat{e}}_{\text{z}},$$(6a) $$\begin{array}{c}{\text{Ω}}_{\text{ss}}\text{}={\text{Ω}}_{\text{ss}}(\frac{\cos {\delta }_{\odot }}{2}(\left(3\text{\hspace{0.17em}}\cos \left(2l^{\prime }-{\alpha }_{\odot }\right)+\cos {\alpha }_{\odot }\right){\widehat{e}}_x\\ +\left(3\text{\hspace{0.17em}Sin}\left(2l^{\prime }-{\alpha }_{\odot }\right)+\sin {\alpha }_{\odot }\right){\widehat{e}}_y)-\sin {\delta }_{\odot }{\widehat{e}}_z),\end{array}$$(6b)

with $${\text{Ω}}_{\text{so}}=\frac{3{(G{M}_{\odot })}^{3/2}}{2c^2{a}^{5/2}{\left(1-e^2\right)}^{5/2}},$$(7a) $${\text{Ω}}_{\text{ss}}=\frac{G{S}_{\odot }}{c^2{a}^3},$$(7b)

where a is the semi-major axis, e is the eccentricity, and f is the true anomaly for Mars. We neglected e in Ω_ss, as the Lense-Thirring effect obtained for a circular orbit is already very small (between three and four orders of magnitude smaller than the geodetic effect; for more, see below).

We now write the components of ê_x, ê_y, and ê_z in the coordinates of a rotating frame attached to Mars (ê_X, ê_Y, ê_Z) and oriented with the Euler angles (ψ, ε, ϕ), with details given in Fig. 1: $$\begin{array}{l}{\widehat{e}}_x=\left(\cos \psi \cos \varphi -\sin \psi \cos \epsilon \sin \varphi \right){\widehat{e}}_X\hfill \\ -\left(\sin \psi \cos \epsilon \cos \varphi +\cos \psi \sin \varphi \right){\widehat{e}}_Y+\sin \psi \sin \epsilon {\widehat{e}}_Z,\hfill \end{array}$$(8a) $$\begin{array}{l}{\widehat{e}}_y=\left(\cos \psi \cos \epsilon \sin \varphi +\sin \psi \cos \varphi \right){\widehat{e}}_X\hfill \\ +\left(\cos \psi \cos \epsilon \cos \varphi -\sin \psi \sin \varphi \right){\widehat{e}}_Y-\cos \psi \sin \epsilon {\widehat{e}}_Z,\hfill \end{array}$$(8b) $${\widehat{e}}_z=\sin \epsilon \sin \varphi {\widehat{e}}_X+\sin \epsilon \cos \varphi {\widehat{e}}_Y+\cos \epsilon {\widehat{e}}_Z.$$(8c)

The precessional angular velocities Ω_so and Ω_ss can also be written as function of the Euler angles such as: $$\begin{array}{l}\Omega =\left(\dot{\epsilon }\cos \varphi +\dot{\psi }\sin \epsilon \sin \varphi \right){\widehat{e}}_X-\left(\dot{\epsilon }\sin \varphi -\dot{\psi }\sin \epsilon \cos \varphi \right){\widehat{e}}_Y\hfill \\ +\left(\dot{\varphi }+\dot{\psi }\cos \epsilon \right){\widehat{e}}_Z,\hfill \end{array}$$(9)

where a "dot" denotes a differentiation with respect to time.

By merging Eq. (9) with Eq. (6a), we obtain $\dot{\epsilon }=\dot{\varphi }=0$, so that only the longitude angle ψ is affected by the geodetic effect $\left(\dot{\psi }={\text{Ω}}_{\text{so}}^z\ne 0\right)$. We proceed to a change of variable (from t to f): $$\text{d}{\psi }_{\text{so}}\left(f\right)={\text{Ω}}_{\text{so}}{\left(1+e\cos f\right)}^3\left(\frac{\text{d}t}{\text{d}f}\right)\text{d}f,$$(10)

with dt/df being given by p^3/2(GM_⊙)^−1/2(1 + e cos f)⁻² for a Keplerian motion, where p = a(1 − e²) the semi-latus rectum of Mars's orbit. After substituting the expression of dt/df into the right-hand side of Eq. (10), the integration is immediate (see also Eq. (3) of Fukushima 1991) and leads to $${\psi }_{\text{so}}\left(f\right)=\frac{3}{2\left(1-e^2\right)}{\left(\frac{na}{c}\right)}^2\left(f+e\sin f\right).$$(11)

After using the equation of the center (see e.g., Murray & Dermott 2000) to express the true anomaly in term of the mean anomaly, l′, the expression for the spin axis longitude of Mars is given by: $$\begin{array}{l}{\psi }_{\text{so}}\left(l^{\prime }\right)=\frac{3}{2\left(1-e^2\right)}{\left(\frac{na}{c}\right)}^2\hfill \\ (l^{\prime }+{\displaystyle \sum _{k=1}^{+\infty }(e\sqrt{1-e^2}\left(j_{k-1}\left(ke\right)-j_{k+1}\left(ke\right)+\frac{2J_k\left(ke\right)}{ke\sqrt{1-e^2}}\right)}\hfill \\ +2{\displaystyle \sum _{m=1}^{+\infty }\frac{{\left(1-\sqrt{1-e^2}\right)}^m}{k{e}^m}}\left(J_{k-m}\left(ke\right)+J_{k+m}\left(ke\right)\right)\sin \left(k{l}^{\prime }\right)),\hfill \end{array}$$(12)

where the J_k(x) are the Bessel functions of first kind with k-index, and n is the Mars's mean motion which is given by Kepler third law of motion. The term which is directly proportional to Mars's mean anomaly describes a precession in longitude at a steady rate, while the other periodic terms represent the nutations in longitude. After making use of the numerical values given in Table 1, we find the following estimate reported here at third-order in eccentricity: $$\begin{array}{ll}{\psi }_{\text{so}}\left(t\right)=\hfill & 6.754\text{\hspace{0.17em}mas\hspace{0.17em}}{\text{yr}}^{-1}\times t+0.565\text{\hspace{0.17em}mas\hspace{0.17em}}\sin l^{\prime }\hfill \\ \hfill & +0.039\text{\hspace{0.17em}mas\hspace{0.17em}}\sin 2l^{\prime }+0.004\text{\hspace{0.17em}mas\hspace{0.17em}}\sin 3l^{\prime }\hfill \end{array}.$$(13)

Since the geodetic precession and nutations are small, the toy model based on the assumption of an elliptic Keplerian orbit is accurate enough for our purposes. The precession and annual terms are above the 0.1 mas threshold and must be included in a model for the longitude angle, ψ, as done, for instance, in Baland et al. (2020). The geodetic precession term is larger than the uncertainty on the determination of the precession rate (−7598.3 ± 2.1 mas yr⁻¹, Le Maistre et al., in prep.) and is needed to avoid an error of about 0.1% in the determination of the polar moment of inertia. The geodetic annual term, with its amplitude of 0.6 mas, does not depend on the properties of Mars's interior and has to be removed from any determination of the annual nutation term before any interpretation in terms the of core radius for instance.

By merging Eq. (9) with Eq. (6b), we obtain the equation of motion for the Euler angles related to the Lense-Thirring effect: $$\begin{array}{ll}{\dot{\psi }}_{\text{ss}}=\hfill & \frac{1}{2}{\text{Ω}}_{\text{ss}}(-2\sin {\delta }_{\odot }\hfill \\ \hfill & +\cos {\delta }_{\odot }\cot \epsilon \left(\sin \left({\alpha }_{\odot }-\psi \right)+3\sin \left(2l^{\prime }-{\alpha }_{\odot }-\psi \right)\right)),\hfill \end{array}$$(14a) $${\dot{\epsilon }}_{\text{ss}}=\frac{1}{2}{\text{Ω}}_{\text{ss}}\cos {\delta }_{\odot }\left(\cos \left({\alpha }_{\odot }-\psi \right)+3\cos \left(2l^{\prime }-{\alpha }_{\odot }-\psi \right)\right),$$(14b) $${\dot{\varphi }}_{\text{ss}}=-\frac{1}{2}{\text{Ω}}_{\text{ss}}\frac{\cos {\delta }_{\odot }}{\sin \epsilon }\left(\sin \left({\alpha }_{\odot }-\psi \right)+3\sin \left(2l^{\prime }-{\alpha }_{\odot }-\psi \right)\right).$$(14c)

After integration (we considered the angles ψ and ε as constant in the right-hand sides of the equations, and therefore denote them with a subscript "0" in the following), the solution for each angle will be the sum of a linear term and of a periodic term at the semi-annual period. After making use of numerical values given in Table 1, we find the following estimate: $$\begin{array}{ll}{\psi }_{\text{ss}}\left(t\right)\hfill & =-0.0779\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\sin {\delta }_{\odot }\hfill \\ \hfill & +0.0390\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\cot {\epsilon }_0\cos {\delta }_{\odot }\sin \left({\alpha }_{\odot }-{\psi }_0\right)\hfill \\ \hfill & -0.0175\text{\hspace{0.17em}μas\hspace{0.17em}}\cos \left(2l^{\prime }-{\alpha }_{\odot }-{\psi }_0\right)\cot {\epsilon }_0\cos {\delta }_{\odot }\hfill \\ \hfill & =-0.0857\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\hfill \\ \hfill & -0.0037\text{\hspace{0.17em}μas\hspace{0.17em}}\cos \left(2l^{\prime }-339°.852\right),\hfill \end{array}$$(15a) $$\begin{array}{ll}{\epsilon }_{\text{ss}}\left(t\right)\hfill & =0.0390\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\cos {\delta }_{\odot }\cos \left({\alpha }_{\odot }-{\psi }_0\right)\hfill \\ \hfill & +0.0175\text{\hspace{0.17em}μas\hspace{0.17em}sin}\left(2l^{\prime }-{\alpha }_{\odot }-{\psi }_0\right)\cos {\delta }_{\odot },\hfill \\ \hfill & =-0.0001\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\hfill \\ \hfill & +0.0017\text{\hspace{0.17em}μas\hspace{0.17em}sin}\left(2l^{\prime }-339°.852\right),\hfill \end{array}$$(15b) $$\begin{array}{ll}{\varphi }_{\text{ss}}\left(t\right)\hfill & =-0.0390\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\frac{\cos {\delta }_{\odot }}{\sin {\epsilon }_0}\sin \left({\alpha }_{\odot }-{\psi }_0\right)\hfill \\ \hfill & +0.0175\text{\hspace{0.17em}μas\hspace{0.17em}sin}\left(2l^{\prime }-{\alpha }_{\odot }-{\psi }_0\right)\frac{\cos {\delta }_{\odot }}{\sin {\epsilon }_0},\hfill \\ \hfill & =-0.0090\text{\hspace{0.17em}μas\hspace{0.17em}}{\text{yr}}^{-1}\times t\hfill \\ \hfill & +0.0040\text{\hspace{0.17em}μas\hspace{0.17em}cos}\left(2l^{\prime }-339°.852\right).\hfill \end{array}$$(15c)

We actually see that the spin-spin contribution is well below the 0.1 mas precision and can therefore be safely neglected, as compared to the spin-orbit contribution.

While considering the spin of Mars as an accelerated gyroscope, we also consider the Thomas precession Ω_ΤΡ contribution within the total relativistic angular precessional velocity Ω (see e.g., Eq. (25) of Soffel et al. 2003). As shown in Eq. (5c), this term is proportional to the acceleration that characterizes the deviation of the actual worldline of the planet from a geodesic, which comes mainly from the coupling of higher order multipole moments of Mars to the external tidal gravitational fields. The Thomas precession scales as Ω_ΤΡ = ‖Ω_TP‖ ∝ Q‖v‖/c², with Q = ‖Q‖ ∝ J₂R²GM_p/‖x − x_p‖⁴; here, J₂ and R are respectively the quadrupole moment and the equatorial radius of Mars, while M_p and x_p are the mass and the barycentric position of the external body, respectively. For the Earth, the dominant contribution to Thomas precession comes from the Moon. Here, for Mars, Phobos plays the dominant role. The absolute value of Q due to the action of Phobos is on the order of 10⁻¹² m s⁻², meaning that Ω_ΤΡ ~ 10⁻⁶ μas yr⁻¹ (against 10⁵ μas yr⁻¹ and 10⁻¹ μas yr⁻¹ for Ω_so and Ω_ss, respectively). The contribution of Thomas precession to the variations in Euler angles is therefore smaller than the already negligible contribution of the Lense-Thirring precession.

2.2 Time coordinate transformation and impact on Mars rotation modeling

We denote, using τ, the local time (or proper time) related to the central body (i.e., Mars) and, using t, the barycentric dynamical time (TDB), which is related to the BCRS and used to analyze the data at Mars. The local model of rotation describes the rotation of the body in terms of local physics such as the atmosphere and surface dynamics. This model is typically expressed as a function of the local time, τ. For example, for a uniformly rotating body, we have $\varphi ={\varphi }_0+{\dot{\varphi }}_0\tau$. As a consequence, when expressed in BCRS, the rotation model will be impacted by the τ − t time transformation. In Sect. 2.2.1, we present the theory related to the time transformation, while in Sect. 2.2.2, we show how this time transformation impacts the modeling of Mars's rotation when expressed in BCRS.

2.2.1 Time coordinate transformation

The t of TDB is a rescaled version of barycentric coordinate time (TCB). The link between TDB and TCB is defined by the recommendation B3 of IAU2006, which reads: $$\frac{\text{d}t}{d\text{TCB}}=1-L_{\text{B}},$$(16)

where L_B = 1.550519768 × 10⁻⁸; as in Eq. (10.3) from Petit et al. (2010).

The proper time related to a local frame co-moving with Mars (the Martian equivalent to geocentric coordinate time, or TCG) is denoted here as τ. For a small orbital velocity (υ ≪ c) and a weak gravitational field (r ≫ GM/c²), the relationship between TDB and Mars's proper time τ at first post-Newtonian order is given by (see, e.g., Soffel et al. 2003; Petit et al. 2010) $$\begin{array}{l}\frac{\text{d}\tau }{\text{d}t}-1={\left[\frac{\text{d}\tau }{\text{d}t}\right]}_{\text{GR}}=\frac{L_{\text{B}}}{1-L_{\text{B}}}-\frac{1}{c^2}\frac{1}{1-L_{\text{B}}}\left(\frac{{\upsilon }^2}{2}+\omega \right)\hfill \\ +\frac{1}{c^4}\frac{1}{1-L_{\text{B}}}\left(-\frac{{\upsilon }^4}{8}-\frac{3}{2}{\upsilon }^2{\omega }^2+4v\cdot w+\frac{1}{2}{\omega }^2+\Delta \right),\hfill \end{array}$$(17)

where v is Mars's barycentric velocity. The potential, w, is the Newtonian potential at the location of Mars: $$w=\left({\displaystyle \sum _{\text{A}\ne \text{Mars}}\frac{G{M}_A}{r_A}}\right)-\frac{3J_2^{\odot }}{2}\frac{G{M}_{\odot }}{r_{\odot }}{\left(\frac{R_{\odot }}{r_{\odot }}\right)}^2\left({\left(\frac{{\widehat{e}}_{\odot }\cdot \left(x-x_{\odot }\right)}{r_{\odot }}\right)}^2-\frac{1}{3}\right),$$(18)

where r_A = ‖x − x_A‖ with x the barycentric position of Mars and x_A is the position of the body A, and where the sum includes the Sun and all planets, R_⊙ is the Sun's equatorial radius, ê_⊙ is the unit vector defining the Sun's spin axis, and $J_2^{\odot }$ is the Sun's quadrupolar moment of the mass distribution. The norm of the barycentric distance between Mars and the Sun is denoted r_⊙ = ‖x − x_⊙‖. In addition, w is the vector potential defined by: $$w={\displaystyle \sum _{\text{A}\ne \text{Mars}}\frac{G{M}_{\text{A}}}{r_A}{v}_{\text{A}},}$$(19)

and Δ is defined by $$\begin{array}{l}\Delta ={\displaystyle \sum _{\text{A}\ne \text{Mars}}\frac{G{M}_{\text{A}}}{r_A}(-2{\upsilon }_{\text{A}}^2}\hfill \\ +{\displaystyle \sum _{\text{B}\ne \text{A}}\frac{G{M}_{\text{B}}}{R_{\text{AB}}}+\frac{1}{2}{\left(\frac{{}_{v_{\text{A}}\cdot \left(x-x_{\text{A}}\right)}}{r_{\text{A}}}\right)}^2}+\frac{1}{2}{a}_{\text{A}}\cdot \left(x-x_{\text{A}}\right)),\hfill \end{array}$$(20)

where a_A is the Newtonian point-mass acceleration of body A, namely, $$a_{\text{A}}={\displaystyle \sum _{\text{B}\ne \text{A}}\frac{G{M}_{\text{B}}}{\left|\right|x_{\text{B}}-x_{\text{A}}|{|}^3}\left(x_{\text{B}}-x_{\text{A}}\right).}$$(21)

We can write the result from the integration of Eq. (17), the relationship between TDB and Mars's propertime, τ, as $$\tau =t+\left[\tau -t\right]\text{GR,}$$(22)

where [τ−t]_GR includes the various relativistic corrections. Typically the relative amplitude of [τ − t]_GR is of the order of GM_⊙/c²a ~ 10⁻⁸, where a is Mars's semi-major axis. In the following, we show how these relativistic corrections impact the data analysis related to the rotation of Mars when analyzed using TDB.

2.2.2 Impact on rotation modeling for Mars

We now explore how the time transformation developed above can impact Mars rotation modeling in BCRS. We first consider a rotation model expressed in Mars's local reference frame. Such a model includes a uniform rotation and periodic terms, that is, $${\varphi }_{\text{local}}\left(\tau \right)={\varphi }_0+{\dot{\varphi }}_{{}_0}^{\text{local}}\tau -{\psi }_{\text{nut}}\cos {\epsilon }_0+{\displaystyle \sum _{j=1}^4\left({\varphi }_{cj}\cos j{l}^{\prime }+{\varphi }_{sj}\sin j{l}^{\prime }\right).}$$(23)

On the other hand, one needs a similar modeling expressed in terms of TDB in order to perform the data analysis. Using Eq. (22), this modeling is given as: $$\begin{array}{l}\varphi \left(t\right)={\varphi }_{\text{local}}\left(\tau \left(t\right)\right)={\varphi }_0+{\dot{\varphi }}_0^{\text{local}}t-{\psi }_{\text{nut}}\cos {\epsilon }_0\hfill \\ +{\displaystyle \sum _{j=1}^4\left({\varphi }_{cj}\cos j{l}^{\prime }+{\varphi }_{sj}\sin j{l}^{\prime }\right)+{\left[\varphi \right]}_{\text{GR}}\left(t\right),}\hfill \end{array}$$(24)

where $${\left[\varphi \right]}_{\text{GR}}\left(t\right)={\dot{\varphi }}_0^{\text{local}}{\left[\tau -t\right]}_{\text{GR}}.$$(25)

When analyzing data from orbiting spacecraft or surface landers, the timescale used in the data reduction is TDB, so that the model for Mars rotation needs to include the relativistic contributions, as presented above. In the following section, we study, both numerically and analytically, various contributions that impact the [τ − t]_GR relationship and their impact on Mars's rotation angle ϕ. We go on to propose a new modeling process that improves upon the one currently used in various analyses.

Additional relativistic corrections in longitude and obliquity can be computed similarly as for the rotation angle, by replacing the rotation rate in Eq. (25) by the precession rate in longitude and in obliquity. However, as those rates (−7598.3 mas yr⁻¹ and −7.9 mas yr⁻¹, respectively, see Le Maistre et al., in prep.) are very small compared to the rotation rate, the associated relativistic corrections are negligible.

3 Impact of the time transformation on Mars rotation modeling

In this section, we study the various contributions to the [τ−t]_GR relationship that arise from the integration of Eq. (17) and their corresponding impacts on the Mars BCRS rotation model through Eq. (25). In order to cross-check and validate our results, we developed three independent approaches.

The first approach consists of developing a simple toy model to analytically integrate Eq. (17). This provides a good physical intuition with regard to the various terms obtained and the process is relatively pedagogical. The second approach consists of numerically integrating Eq. (17) using the DE440 planetary ephemerides (Park et al. 2021) provided by the NAIF-SPICE software (Acton 1996; Acton et al. 2018). In a subsequent step, analytical series consisting in various harmonic terms are fitted to the result of the numerical integration. This procedure is similar to the one developed to produce the Ephemeris Time (ET) for Earth. For details, we refer to Fukushima (1995, 2010), Irwin & Fukushima (1999), and Harada & Fukushima (2003).

The third approach we developed consists of using a series for the barycentric position and distance of planets to obtain a subsequent series for the relativistic part of the rotation angle. We use the analytical planetary theory VSOP87 (Bretagnon & Francou 1988), derived from the DE200 planetary ephemerides (Standish 1982). Using the VSOP87 series, we analytically integrate Eq. (17) and identify the harmonics with the largest contribution to the [τ−t]_GR relationship. We purposely used VSOP87 in this instance, instead of the more recent versions of VSOP theories that are not suited for our purpose. VSOP2000 (Moisson & Bretagnon 2001) provides a series only for the heliocentric positions (not the barycentric positions), whereas VSOP2013 (Simon et al. 2013), derived from the INPOP planetary ephemerides (Fienga et al. 2011; Bernus et al. 2019, 2022), is based on Tchebychev polynomials (and not on a series). We go on to show that the error introduced by using the older (and therefore less accurate) VSOP87 theory is negligible for our purpose, as the solutions of the second and third approaches are consistent to the 0.1 mas level.

The advantages to considering the last two approaches are rooted in the fact that we can cross-check our results and estimate the uncertainties in our derived modeling coming from numerical integration, differences between the DE and INPOP ephemerides, and so on.

In the remainder of this section, we consider the various contributions from Eq. (17): the 1/c² contribution coming from the two-body problem, contribution related to the motion of the Sun with respect to the Solar System barycenter (SSB), contribution from the L_B constant, direct contribution from other planets, and higher order contributions (i.e., the 1/c⁴ terms and the Sun's $J_2^{\odot }$).

3.1 Simple analytical solution (toy model)

3.1.1 The 1/c² contribution coming from the Sun considering a Keplerian motion

We first consider the main contribution to the relation between τ and TDB, which is the 1/c² contribution from the Sun in a two-body (or Keplerian) problem.

The evolution of proper time with respect to coordinate time is given by: $$\frac{\text{d}\tau }{\text{d}t}-1={\left[\frac{\text{d}\tau }{\text{dt}}\right]}_{\text{2body}}=\frac{L_{\text{B}}}{1-L_{\text{B}}}-\frac{1}{1-L_{\text{B}}}\frac{{\upsilon }_{M\odot }^2}{2c^2}-\frac{1}{1-L_{\text{B}}}\frac{G{M}_{\odot }}{r_{M\odot }{c}^2},$$(26)

where r_M⊙ is the distance between the Sun and Mars and v_M⊙ is the norm of their relative velocity. This expression can be integrated exactly assuming a perfect Keplerian motion, as in, for instance, Moyer (1981): $$\begin{array}{ll}{\left[\tau -t\right]}_{\text{2body}}\hfill & =\text{cst}+\frac{L_{\text{B}}}{1-L_{\text{B}}}t-\frac{1}{1-L_{\text{B}}}\frac{n{a}^2}{2c^2}\left(4E-nt\right),\hfill \\ \hfill & =\text{cst}+\frac{L_{\text{B}}}{1-L_{\text{B}}}t-\frac{1}{1-L_{\text{B}}}\frac{n{a}^2}{c^2}\left(2e\sin E+\frac{3}{2}{l}^{\prime }\right),\hfill \end{array}$$(27)

where n is the mean motion, a the semi-major axis, E is the eccentric anomaly, and l′ is the mean anomaly. A low eccentricity expansion leads to: $$\begin{array}{ll}{\left[\tau -t\right]}_{\text{2body}}=\hfill & \text{cst}+\frac{1}{1-L_{\text{B}}}\left(L_{\text{B}}-\frac{3n^2{a}^2}{2c^2}\right)t\hfill \\ \hfill & -\frac{n{a}^2}{c^2\left(1-L_{\text{B}}\right)}\left(2e-\frac{e^3}{4}\right)\sin l^{\prime }\hfill \\ \hfill & -\frac{n{a}^2}{c^2\left(1-L_{\text{B}}\right)}\left(e^2-\frac{e^4}{3}\right)\sin 2l^{\prime }\hfill \\ \hfill & -\frac{3n{a}^2}{4c^2\left(1-L_{\text{B}}\right)}{e}^3\sin 3l^{\prime }\hfill \\ \hfill & -\frac{2n{a}^2}{3c^2\left(1-L_{\text{B}}\right)}{e}^4\sin 4l^{\prime }+\mathrm{...},\hfill \end{array}$$(28)

which includes a linear drift and oscillations at frequencies multiple of the orbital frequency. The linear drift includes a contribution from the rescaling between TCB and TDB (i.e., the contribution from L_B) of 1.55 × 10⁻⁸ and a contribution of −9.72 × 10⁻⁹ from the Sun (using Mars's orbital parameters from Table 1). The total linear drift coefficient is 5.79 × 10⁻⁹. The amplitudes of the harmonics terms are: −11.419 ms for the term at orbital period, −532.3 μs for the term at twice the orbital period, −37.4 μs for the term at three times the orbital period, and −3.1 μs for the term at four times the orbital period. We note that the 1/(1 − L_B) coefficient impacts these amplitude only at the relative level of 10⁻⁸.

The impact on Mars's rotation is obtained from Eq. (24). First of all, it is important to notice that the ${\dot{\varphi }}_0$ estimated from a data analysis performed in BCRS is actually: $${\dot{\varphi }}_0={\dot{\varphi }}_0^{\text{local}}\left(1+\frac{1}{1-L_{\text{B}}}\left(L_{\text{B}}-\frac{3n^2{a}^2}{2c^2}\right)\right)={\dot{\varphi }}_0^{\text{local}}\left(1+5.79\times {10}^{-9}\right),$$(29)

where ${\dot{\varphi }}_0^{\text{local}}$ is the proper Mars rotation rate. Using the measured value of ${\dot{\varphi }}_0$ from Table 1, we can find: $${\dot{\varphi }}_0^{\text{local}}=350°.891983308{\text{day}}^{-1},$$(30)

such that ${\varphi }_{\text{GR}}=5.79\times {10}^{-9}{\dot{\varphi }}_0^{\text{local}}$ is 2°.03 × 10⁻⁶ day⁻¹ (or 7.3117 mas day⁻¹). This quantity is two orders of magnitude larger than the current uncertainty in the Mars rotation rate estimate and thus must be removed for any geophysical interpretation of the latter (${\dot{\varphi }}_0$ in Eq. (2) must be replaced by ${\dot{\varphi }}_0^{\text{local}}$).

The 2-body contribution from the transformation between τ and t to the Mars rotation model in BCRS, including the linear term and the four largest periodic terms, is $$\begin{array}{lll}{\left[\varphi \right]}_{2\text{body}}\hfill & =\hfill & {\dot{\varphi }}_0^{\text{local}}{\left[\tau -t\right]}_{\text{GR}}={\dot{\varphi }}_{\text{GR}}t+{\displaystyle \sum _j{\varphi }_{rj}\sin j{l}^{\prime },}\hfill \\ \hfill & =\hfill & 2°.03\times {10}^{-6}{\text{day}}^{-1}\times t\hfill \\ \hfill & \hfill & -166.950\text{\hspace{0.17em}mas\hspace{0.17em}}\sin l^{\prime }-7.782\text{\hspace{0.17em}mas\hspace{0.17em}}\sin 2l^{\prime }\hfill \\ \hfill & \hfill & -0.547\text{\hspace{0.17em}mas\hspace{0.17em}}\sin 3l^{\prime }-0.045\text{\hspace{0.17em}mas\hspace{0.17em}}\sin 4l^{\prime }.\hfill \end{array}$$(31)

These values will be refined in the further subsections considering a more accurate modeling of Mars's trajectory. Nevertheless, the Keplerian modeling presented here is sufficient to get an estimate of the order of magnitude of the impact of the time transformation on the rotation of Mars.

3.1.2 Contribution related to the motion of the Sun with respect to the Solar System barycenter

The two-body problem calculation performed in the previous section considers one test mass orbiting one massive body, assuming that the coordinate time TDB is the one related to a coordinate system where the massive body is at rest at the origin. For BCRS, this is not the case: the Sun is not at rest and not located at the origin of the coordinate system that is defined as the SSB. Therefore, when we are interested in computing the evolution of Mars proper time with respect to TDB, we should use: $$\frac{d\tau }{dt}-1=\frac{L_{\text{B}}}{1-L_{\text{B}}}-\frac{1}{1-L_{\text{B}}}\frac{{\upsilon }_{\text{M}}^2}{2c^2}-\frac{1}{1-L_{\text{B}}}\frac{G{M}_{\odot }}{r_{\text{M}\odot }{c}^2},$$(32)

where v_M is Mars's velocity with respect to the SSB and r_M⊙ is the distance between Mars and the Sun. The only difference with respect to Eq. (26) relies on the velocity expressed relative to the SSB and not to the Sun. A simple calculation using v_M = v_M⊙ + v_⊙, where v_M⊙ is the velocity of Mars with respect to the Sun and v_⊙ is the Sun velocity with respect to the SSB, shows that there is an additional contribution to the evolution of Mars proper time due to the velocity of the Sun with respect to the SSB. This additional contribution is given by: $${\left[\frac{d\tau }{dt}\right]}_{\text{SSB}}=-\frac{1}{1-L_{\text{B}}}\frac{{\upsilon }_{\odot }\cdot v_{\text{M}\odot }}{c^2}-\frac{1}{1-L_{\text{B}}}\frac{{\upsilon }_{\odot }^2}{2c^2}\approx -\frac{v_{\odot }\cdot v_{\text{M}\odot }}{c^2}.$$(33)

To first order, the motion of the Sun with respect to the SSB is due to its gravitational interaction with Jupiter. If we take a simple toy model and consider that the motion of the Sun with respect to the SSB is due to Jupiter only and assuming Jupiter's orbit to be circular around the Sun, then the Sun's velocity is given by ${\upsilon }_{\odot }\approx a_{\text{J}}{n}_{\text{J}}\frac{M_{\text{J}}}{M_{\odot }}$ (sin n_Jt, − cos n_Jt, 0), where a_J is Jupiter's semi-major axis and n_J its mean motion. To a first approximation, we can also consider the orbital motion of Mars to be circular; then, Eq. (33) can be integrated analytically: $${\left[\tau -t\right]}_{\text{SSB}}\approx \frac{a{a}_{\text{J}}}{c^2}\frac{n{n}_{\text{J}}}{n-n_{\text{J}}}\frac{M_{\text{J}}}{M_{\odot }}\sin \left(\left(n-n_{\text{J}}\right)t+\delta \right),$$(34)

where a is Mars's semi-major axis, n its mean motion, and δ the phase difference between the two planets, assuming co-planar motion. The velocity of the Sun with respect to the SSB therefore induces an additional modulation to the transformation from τ to the TDB. The period of this modulation is the Mars-Jupiter synodic orbital period, namely, 2.235 yr, and its amplitude is 37.59 μs.

Using Eq. (25), this modulation impacts the BCRS Mars rotation modeling and induces a modulation of amplitude of 0.55 mas at the Mars-Jupiter synodic orbital period. A similar calculation considering Saturn leads to a harmonic term of period of 2 yr with an amplitude of 0.11 mas. The other planets induce periodic terms with amplitudes ≲0.01 mas. This can be seen as an indirect effect of the other planets of the Solar System on Mars propertime since it comes from the impact of other planets on the SSB velocity. The direct effect will be computed below.

3.1.3 Direct contribution from other planets

As can be noticed from Eq. (17), the gravitational potential from the other planets will also impact the evolution of Mars proper time. To first order, the impact from the planets gravitational potential is governed by $${\left[\frac{d\tau }{dt}\right]}_{\text{P}}=-\frac{1}{1-L_{\text{B}}}\frac{G{M}_{\text{P}}}{c^2{r}_{\text{M}P}},$$(35)

where r_MP = ‖x_M − x_P‖ is the distance between Mars and a planet P.

A simple toy model considering both Mars and the planet to be orbiting on coplanar circular orbits shows that, to first order, the integration of the previous equation leads to a linear drift whose linear coefficients is given by $-G{M}_{\text{P}}/\left({\left(a^2+a_{\text{P}}^2\right)}^{1/2}{c}^2\right)$ and to an harmonic signal at the planet-Mars synodic orbital period and of amplitude of $G{M}_{\text{P}}\frac{a{a}_{\text{P}}}{{\left(a^2+a_{\text{P}}^2\right)}^{3/2}}/\left(c^2\left(n-n_{\text{P}}\right)\right)$, where a_P and n_P are the semi-major axis and mean motion of the planet P. This calculation is valid only to first order in $a{a}_{\text{P}}/\left(a^2+a_{\text{P}}^2\right)$ and neglecting L_B. Other harmonics can be identified at higher orders and for non-zero eccentricities, in particular an oscillation at the orbital period of the planet (see Sect. 3.3).

For each planet P, these contributions to the τ to TDB transformation will impact Mars's BCRS rotation modeling through Eq. (25), that is to say, it will produce one term with linear rate ${\dot{\varphi }}_{\text{GR}}$ and one harmonic synodic term. For Jupiter, ${\dot{\varphi }}_{\text{GR}}=-6°\times {10}^{-10}{\text{day}}^{-1}$ (or −0.00220 mas day⁻¹). For Saturn, ${\dot{\varphi }}_{\text{GR}}=-1°\times {10}^{-10}{\text{day}}^{-1}$ (or −0.00037 mas day⁻¹). The synodic terms associated to Jupiter and Saturn have amplitudes of 0.077 mas and 0.007 mas, respectively. The other planets induce synodic terms with amplitudes ≲0.001 mas. Although the direct synodic terms associated to Jupiter and Saturn have small amplitudes, they cannot be neglected as they combine with the indirect synodic terms obtained in the previous subsection. The total synodic terms related to Jupiter and Saturn have 0.47 mas and 0.10 mas of amplitude, respectively, as the direct and indirect terms are out of phase to each other. The sum of the direct and indirect effect of the planets will be refined in the further subsections considering a more accurate modeling of each planet’s trajectory.

Similarly as for the effect of the planets, it is possible to build a toy model for the direct effect of Phobos and Deimos, and of Ceres, the largest body of the asteroid belt. However, given their small mass, the associated relativistic corrections can be neglected.

3.2 Numerical solution using the DE planetary ephemerides

In this section, we present the result of a numerical integration of Eq. (17) using the DE440 planetary ephemerides (Park et al. 2021) provided by the NAIF-SPICE software (Acton 1996; Acton et al. 2018). The integration was performed starting from J2000 and lasting 30 yr backward and forward. In a second step, analytical series consisting of various harmonic terms were fitted to the result of the numerical integration. This procedure is similar to the one developed to produce the Time Ephemeris for Earth (see Fukushima 1995, 2010; Irwin & Fukushima 1999; Harada & Fukushima 2003).

3.2.1 First-order contributions

In this section, we consider the leading contributions from Eq. (17), namely, the 1/c² contributions from the Sun, the various planets, and from L_B. This integration therefore includes all the effects presented in the previous section. We numerically integrate Eq. (17) and transform the evolution of τ − t into an estimate of the evolution of ϕ through Eq. (25). We then fit the expression $${\left[\varphi \right]}_{\text{GR}}\left(t\right)={\dot{\varphi }}_{\text{GR}}t+{\displaystyle \sum _j\left(C_j\cos f_{j}t+S_j\sin f_{j}t\right)}$$(36)

to the numerically integrated evolution of [ϕ]_GR(t). The values of the various coefficients ${\dot{\varphi }}_{\text{GR}}$, C_j, and S_j are obtained using a standard linear least-squares fit. Motivated by the toy model presented in the previous section, the angular frequencies, f_j, included in the fit are chosen as linear combinations of the planets mean motion. We identify the relevant frequencies by (iteratively) searching for the largest peaks in the Fourier transform of the numerically integrated [ϕ]_GR(t) time series. The fitted coefficients are then transform to obtain the following expression: $${\left[\varphi \right]}_{\text{GR}}\left(t\right)={\dot{\varphi }}_{\text{GR}}t+{\displaystyle \sum _j{A}_j\sin \left(f_{j}t+{\phi }_j\right).}$$(37)

The estimated linear term is ${\dot{\varphi }}_{\text{GR}}=2°.03021\times {10}^{-6}{\text{day}}^{-1}$ (or 7.308758 mas day⁻¹), while the amplitude and phase of the harmonic terms are given in Table 2. We searched for terms with amplitude down to 0.04 mas, covering the amplitude range considered in Eq. (31) for the terms at the harmonics of the Martian orbital period. For these terms, the estimated solution is in good agreement with the solution of the toy model (difference <0.005 mas), but differs significantly with respect to Eq. (21) of Yoder & Standish (1997), recalled in Eq. (3), with up to 9 mas, or 5% in annual amplitude. The differences are mainly likely due to truncation errors in the parameters values used by Yoder & Standish (1997). The estimated terms at the Mars-Jupiter and Mars-Saturn synodic periods are 0.10 mas and 0.005 mas larger than obtained with the toy model, as a result of the assumption of circular planetary orbits therein. We also find one term at the orbital period of Jupiter with an amplitude of about 0.04 mas, related to the direct effect of the planet, and three others terms at different periods with amplitude ranging between 0.06 mas and 0.08 mas, mainly due to the indirect effect of the Earth and of Jupiter of the orbit of Mars.

The residuals between the numerical integration and the fitted harmonics decomposition is presented in Fig. 2 and remains below 0.15 mas. Formal uncertainties are not relevant quantities to characterize the errors of the fit as they do not directly rely on any observations (no data points). Instead, we use another method described in Sect. 3.3.2 to compute the accuracy of the estimated coefficients, which is equal to ~0.01 mas, on the main terms.

Fig. 2 Difference between the [ϕ]GR(t) obtained by numerically integrating Eq. (17), using the DE440 planetary ephemerides and using Eq. (25) and the fitted series from Eq. (37). This curve provides an estimate of the accuracy of the fitted analytical model provided by the coefficients from Table 2.

Fig. 3 Impact of the 1/c4 term from Eq. (17) on the BCRS modeling of Mars’s rotation (top). Impact of the Sun quadrupole moment $\left(J_2^{\odot }\right)$ on the BCRS modeling of Mars’s rotation (bottom).

3.2.2 Higher-order contributions

In this section, we consider the 1/c⁴ contribution appearing in Eq. (17) and the contribution from the Sun’s quadrupole moment, $J_2^{\odot }$. Their impact on the BCRS modeling of Mars’s rotation is presented in Fig. 3 and is on the order of μ as for the 1/c⁴ term and on the order of 10 nas for the $J_2^{\odot }$. Both these contributions can safely be neglected.

3.3 Series solution using VSOP ephemerides

In this section, we present the results of a semi-analytical approach based on series for the barycentric position and distance of planets as provided by the analytical planetary theory VSOP87 (Bretagnon & Francou 1988). We integrate Eq. (17), neglecting the 1/c⁴ and Sun’s $J_2^{\odot }$ contributions, identify the harmonics with the largest contribution to the τ − t relationship, and estimate the evolution of ϕ through Eq. (25). In a second step, we numerically assess the accuracy of the semi-analytical solution.

3.3.1 Semi-analytical solution

In VSOP87, the barycentric Cartesian coordinates (X, Y, Z) and distance r of the planets to the Sun are written as series of the form: $${\displaystyle \sum _j\left(C_j\cos {\phi }_j+S_j\sin {\phi }_j\right),}$$(38a) $$C_j={\displaystyle \sum _{\alpha }{T}^{\alpha }{C}_{\alpha ,j},}$$(38b) $$S_j={\displaystyle \sum _{\alpha }{T}^{\alpha }{S}_{\alpha ,j},}$$(38c)

where C_α,j and S_α,j are amplitudes and T is the time measured in thousands of Julian years from J2000. φ_j are linear combinations of fundamental arguments, including the mean longitudes of Saturn (Sa), Jupiter (Ju), Mars (Ma), and the Earth (Te); we refer to Table 2 of (Bretagnon & Francou 1988): $$Sa\text{ = 0.87401675650 + 213.2990954380 }T\text{,}$$(39a) $$Ju\text{ = 0.59954649739 + 529.6909650946 }T\text{,}$$(39b) $$Ma\text{ = 6.20347611291 + 3340.6124266998 }T\text{,}$$(39c) $$Te\text{ = 1.75347045953 + 6283.0758499914 }T\cdot$$(39d)

The power α is an integer in-between 0 and 5. For α = 0, the series are periodic. For α ≥ 1, the series are pseudo-periodic (Poisson series).

The solution for [ϕ]_GR(t) is firstly written as: $${\left[\varphi \right]}_{\text{GR}}\left(t\right)\text{=}{\dot{\varphi }}_{\text{GR}}\times t+{\displaystyle \sum _j\left(\Delta {\varphi }_j^c\cos {\phi }_j+\Delta {\varphi }_j^s\sin {\phi }_j\right),}$$(40a) $$\Delta {\varphi }_j^{c/s}\text{=}{\displaystyle \sum _{\alpha }{T}^{\alpha }\Delta {\varphi }_{\alpha ,j}^{c/s},}$$(40b)

with $\text{Δ}{\varphi }_j^{c/s}$ the amplitudes of the periodic and Poisson series. In the first place, the fundamental arguments of the series will be the same as the VSOP87 arguments, by construction, because the series for υ², the squared Mars barycentric velocity, is directly obtained as the squared norm of the time derivative of the position vector of Mars (X, Y, Z), while the series for 1/r_⊙ is obtained starting with the VSOP series for r_⊙ and following Eq. (61) of Baland et al. (2020). The distance r_P between Mars and another planet varies greatly with time, and as a result, it is difficult to express as convergent series for 1/r_P starting from the VSOP series for the Cartesian coordinates. We therefore assume that the orbits of Mars and of the other planets are Keplerian and coplanar and use the mean orbital elements of Simon et al. (2013). We adapted the procedure described in Sect. 4.3.3 of Baland et al. (2020) to obtain a series for $1/r_p^5$ to the case of 1/r_P. This can be seen as extension of the toy model presented in Sect. 3.1.3 to higher orders in eccentricities and in $a{a}_p/\left(\right(a^2+a_{\text{p}}^2\left)\right)$ and, as a result, we go on to identify more harmonics and obtain different amplitudes.

Then, for consistency with the usual form of [ϕ]_GR(t), expressed with the mean anomaly, l′, of Mars as the argument (see Eq. (3)), we change the fundamental arguments, using the mean anomalies of the planets $$l_{Sa}=Sa-{\varpi }_{Sa}=5\text{.53304687684 + 213.2002152909 }T\text{,}$$(41a) $$l_{Ju}=Ju-{\varpi }_{Ju}=0\text{.52395267692 + 529.6533496052 }T\text{,}$$(41b) $$l^{\prime }=Ma-{\varpi }_{\text{M}a}=0\text{.3381185455 + 3340.5349512479 }T\text{,}$$(41c) $$l_{Te}=Te-{\varpi }_{Te}=6\text{.24006011944 + 6283.0195517158 }T\text{,}$$(41d)

instead of their mean longitudes. We express [ϕ]_GR(t) correct up to the first order in the rates of the pericenter longitudes ϖ of the planets, creating a second Poisson series, to add to the first Poisson series coming directly from VSOP ephemerides, and which is not affected by the argument change, at first order. Both Poisson series are similar, but with opposite amplitudes, and therefore they almost cancel each other (the sum of the two series is smaller than 0.05 mas on the interval ±30 yr around J2000 and significant Poisson terms were not found in the fit of the numerical solution). As a result, we omitted Poisson series in the following.

Finally, the periodic series in [ϕ]_GR(t) is written in a pure Sine form, convenient for the purposes of application: $${\left[\varphi \right]}_{\text{GR}}\left(t\right)={\dot{\varphi }}_{\text{GR}}t+{\displaystyle \sum _j{\varphi }_j^r\sin \left(f_{j}t+{\phi }_j^0\right),}$$(42)

with ${\varphi }_j^r$ the amplitudes and ${\phi }_j^0$ a phase (different from that of Eq. (40a)); f_j represents the linear combination of the rate of the mean anomalies of Eq. (41).

For the constant rate, we find ${\dot{\varphi }}_{\text{GR}}=7.3088$ mas day⁻¹, in agreement with the sum of the respective contributions from L_B, from the Sun and from each planet, as obtained with the fit of the numerical solution. The periodic terms, with amplitudes down to 0.04 mas, are presented in Table 3. We find the same ten periodic terms as with the fit of the numerical solution of Table 2, but also one additional term, at the orbital period of Saturn (~30 yr) and with an amplitude of 0.043 mas. We did not find this term with the fit because of its long period and of its small amplitude.

Fig. 4 Difference between the solutions for [ϕ]GR(t) obtained using VSOP87 and DE440 ephemerides via numerical integrations (left). Difference between the numerical integration and the semi-analytical solution using VSOP87, with only the 11 terms of Table 3 in blue, with all terms in black (right).

3.3.2 Accuracy of the semi-analytical solution

The difference between the numerical integrations of [ϕ]_GR(t) performed using the recent DE440 and the older VSOP87 ephemerides is presented in the left panel of Fig. 4, and it remains below 0.003 mas, indicating that using the VSOP87 theory should not be a cause of major errors.

The different steps of the computational procedure to obtain the semi-analytical solution of Eq. (40a) introduce residuals smaller than 0.05 mas when all the terms of the periodic series are considered, and smaller than 0.2 mas when only the 11 largest terms of the series are considered (see right panel of Fig. 4). This is in agreement with the residuals of the fitted solution to the numerical integration based DE440 ephemerides, which includes the ten periodic terms of Table 2 (see Fig. 2).

The amplitudes and phases of the fitted (Table 2) and semi-analytical (Table 3) solutions are in good agreement to each other (see also Fig. 5, where the difference remains below 0.06 mas). The differences give a sense of the modeling uncertainties: for instance, 0.01 mas (0.005%) in annual amplitude or 0.01 mas (0.9%) on the Mars-Jupiter synodic term. Both solutions or an averaged solution can be used for the purposes of application.

4 Signatures in the Doppler observable

As shown in the previous sections, the relativistic variations in Mars rotation and orientation mainly affect the angles ψ and ϕ. Using analytic expressions, we characterize in this section the signature of these variations in the Doppler observable of a Martian lander communicating directly with the Earth.

Numerical applications are provided for the specific case of RISE (Rotation and Interior Structure Experiment), the radio-science experiment of the NASA InSight mission (Folkner et al. 2018). The level of the signatures and their temporal behavior are compared to the noise level and non-relativistic signatures, respectively.

We note $\delta \left({\dot{\psi }}_{0}t\right)$ as a variation in the precession, δ(ψ_nut) a variation in the nutation in longitude, and δϕ a variation (linear and/or periodic) in the rotation angle ϕ, excluding the nutation term ψ_nut cos ε₀ (see Eqs. (1) and (2)). These three notations can be used to represent any kind of variations in ψ and ϕ, including the relativistic ones. Considering that #3, δ(ψ_nut), and δϕ are small quantities, their signature in the observable can be written to the first order as (Yseboodt et al. 2017): $$\begin{array}{ll}\Delta q_{\delta \left({\dot{\psi }}_{0}t\right)}\hfill & =\delta \left({\dot{\psi }}_{0}t\right)\text{Ω}R\cos \theta \hfill \\ \hfill & \left(\cos {\delta }_{\text{E}}\cos H_{\text{E}}\cos \epsilon -\sin {\delta }_{\text{E}}\sin \left(H_{\text{E}}+{\alpha }_{\text{E}}\right)\sin \epsilon \right),\hfill \end{array}$$(43a) $$\Delta q_{\delta \left({\psi }_{\text{nut}}\right)}=-\delta \left({\psi }_{\text{nut}}\right)\text{Ω}R\cos \theta \sin {\delta }_{\text{E}}\sin \left(H_{\text{E}}+{\alpha }_{\text{E}}\right)\sin \epsilon ,$$(43b) $$\Delta q_{\delta \varphi }=\delta \varphi \text{\hspace{0.17em}Ω}R\cos \theta \cos H_{\text{E}}\cos {\delta }_{\text{E}},$$(43c)

with θ as the lander latitude, R the radius of Mars, Ω the rotation rate, and δ_E the Earth declination relative to Mars’s equator. Each of these Doppler variations has a diurnal modulation via, H_E, the Earth hour angle seen from Mars. The observable variations defined in Eq. (43) are Doppler shift expressed as the variation of the velocity along the line of sight (LOS). For a round-trip (two-way) radio link, the conversion factor between the LOS velocity and the Doppler observable is 2f_t/c, with f_t as the carrier downlink frequency. In the case of RISE (but also of LaRa, a radio transponder ready to fly to Mars, see Dehant et al. 2020), f_t ≃ 8.4 GHz, which is in the X frequency band.

The signature of [ϕ]_GR(t) (linear and/or periodic terms) is computed using Eq. (43c). The signature of the geodetic precession in longitude is obtained using Eq. (43a) and that of the geodetic nutation in longitude is computed with Eq. (43b). Precession and nutation signatures differ from each other because ψ_nut also affects the angle ϕ, while ${\dot{\psi }}_0$ does not (see Eqs. (1) and (2)).

The periodic variations in ϕ(t) induced by the seasonal atmosphere and surface dynamics (see, e.g., Konopliv et al. 2020) and by the time coordinates transformation altogether result in a maximum angular displacement of the lander of ~670 mas as seen from the center of Mars (1 mas corresponds to a displacement of 1.6 cm at the surface of Mars). For a lander located at the InSight landing site (i.e., Elysium Planitia: 4.5°N, 135.62°E, −2.6 km altitude), such an angular displacement induces a Doppler shift in the RISE measurements of ≲0.56 mm s⁻¹ (see Fig. 2 of Yseboodt et al. 2017 or Table 4). A fourth of this Doppler signal, computed at the RISE tracking data-timing using the RISE frequency, comes from the relativistic periodic variations (see Fig. 6a and extra information in Table 4).

The combination of seasonal variations in Euler angles and of diurnal trend in the hour angle produces symmetrical envelopes with respect to zero in the Doppler signature, as seen in Figs. 6 and 7. Because of the repeatability in the RISE observation timing imposed by its fixed and directional antennas, the data points cover a limited part of these diurnal cycles.

The ~10 mas difference between the periodic terms of our solution for [ϕ]_GR and that first estimated by Yoder & Standish (1997, see Eq. (3) and Sect. 3.2.1) has a Doppler signature lower than 0.007 mm s⁻¹, which is smaller than the RISE noise level (1.1 mHz at 60 s integration time, corresponding to 0.02 mm s⁻¹) and smaller than the liquid core signature (~0.01 mm s⁻¹, Le Maistre et al., in prep). However, a precise solution for the Mars rotation angles (at the level of 1 mas or smaller) is needed to correctly interpret the measured periodic variations in terms of atmospheric constraints.

The signature in the Doppler of the relativistic linear term of 7.3 mas day⁻¹ (see Eq. (46c)) is very large (up to 49 mm s⁻¹), as shown in Fig. 6b. This pronounced signal with a linear increase that is barely visible in the plot occurs because it has a linear dependence on the time of the observations (2018–2022 for InSight) relative to the chosen reference epoch (here, J2000).

The signature of the geodetic nutation in longitude in the Doppler observable (shown in Fig. 7a) is very small (≤ 0.12 µm s⁻¹), while that of the geodetic precession (linear term in Eq. (13)) is two to three orders of magnitude larger (up to 0.03 mm s⁻¹), as shown in Fig. 7b, where the signature of the linear plus periodic geodetic terms is plotted (see Eq. (13)). Similarly as for the relativistic linear term in the rotation angle, the stronger signature of the geodetic linear term in longitude for Mars (6.75 mas yr⁻¹) results from the linear dependency to the time past from the chosen origin of time (i.e., ~150 mas in 2020 for a reference epoch at J2000).

Fig. 5 Difference between the [ϕ]GR(t) of the semi-analytical solution using VSOP87 (with only the 11 terms of Table 3) and the fit of the numerical integration with DE440 ephemerides.

5 Application to other planets

Mars has been extensively explored, but the rotation of other bodies of the Solar System is also subject to investigations. Here, we provide an estimate of the main terms of the relativistic correction in longitude, ψ, and rotation, ϕ, to include in the rotation model of our neighboring planets. For each planet, we first estimated the geodetic precession and nutation in longitude. Then we estimated the linear and periodic changes in rotation, considering (1) the Sun’s contribution with regard to a planet on a Keplerian orbit, (2) the contribution related to the Sun motion with respect to the SSB, and (3) the direct contribution from the other planets.

5.1 Geodetic precession and nutations

For any planet, the spin–orbit angular velocity Ω_so of Eq. (5a) is proportional to x × v and therefore perpendicular to the orbital plane (see also Eq. (6a)). As a result, if a planet moves on a Keplerian orbit, only its longitude angle (defined as the longitude of the spin axis with respect to the orbital plane) is affected by the geodetic precession and nutation, whereas the obliquity and rotation angles ε and φ are unchanged: $${\dot{\psi }}_{\text{so}}={\text{Ω}}_{\text{so}}^z,$$(44a) $${\dot{\epsilon }}_{\text{so}}=0,$$(44b) $${\dot{\varphi }}_{\text{so}}=0.$$(44c)

The geodetic precession and nutation for each planet, as obtained with the toy model of Sect. 2.1, is given in Table 5. The geodetic precession rate increases with decreasing distance to the Sun. The amplitudes of the periodic terms do not strictly follow that rule of thumb because they are relatively more dependent on eccentricity.

For the Earth, we obtain consistent results with Fukushima (1991), as we follow the same approach based on a Keplerian orbit (see also Eq. (27) of Soffel et al. 2003). Even though their solution is presented as an approximation, our results for the precession and annual terms of the eight planets are in a very good agreement with those presented in Table 1 of Eroshkin & Pashkevich (2007), where they fit a solution to a numerical integration based on ephemerides. This is because of two approximations which compensate each other during their computation: (1) they refer the geodetic motion of all planets to the Earth ecliptic of J2000, instead of their respective orbital plane, and (2) they neglect the equatorial components of the angular velocity vector σ = (σ^X, σ^Y‚ σ^Z) expressed in the coordinates of a frame attached to the Earth ecliptic of J2000. For the demonstration, we first write the geodetic variations in “ecliptic Euler angles” (we use the notation * for these angles) as $${\dot{\psi }}_{\text{so}}^{*}={\sigma }^Z-{\sigma }^X\sin {\psi }_0^{*}\cot {\epsilon }_0^{*}+{\sigma }^Y\cos {\psi }_0^{*}\cot {\epsilon }_0^{*},$$(45a) $${\dot{\epsilon }}_{\text{so}}^{*}={\sigma }^X\cos {\psi }_0^{*}+{\sigma }^Y\sin {\psi }_0^{*},$$(45b) $${\dot{\varphi }}_{\text{so}}^{*}=\left({\sigma }^X\sin {\psi }_0^{*}-{\sigma }^Y\cos {\psi }_0^{*}\right)/\sin {\epsilon }_0^{*}.$$(45c)

Then we write σ = R_z(−Ω₀) · R_x(−i₀) · Ω_so, with Ω₀ and i₀ the ascending node longitude and inclination of the planet’s orbit with respect to the ecliptic. Since ${\text{Ω}}_{\text{SO}}^x\simeq 0$ and ${\text{Ω}}_{\text{SO}}^y\simeq 0$ (for a non Keplerian orbit, small variations about zero are possible), ${\sigma }^X\simeq {\text{Ω}}_{SO}^z\sin i_0\sin {\text{Ω}}_0$, ${\sigma }^Y\simeq -{\text{Ω}}_{\text{SO}}^z\sin i_0\cos {\text{Ω}}_0$, and ${\sigma }^Z\simeq {\text{Ω}}_{\text{SO}}^z\cos i_0$.

To obtain ψ_so(t), we integrate ${\text{Ω}}_{\text{SO}}^z$ over time (see Eq. (44a)). For bodies with a small orbital inclination with respect to the Earth’s ecliptic, ${\sigma }^Z\simeq {\text{Ω}}_{\text{SO}}^z$, and by integrating only σ^Z while neglecting σ^X and σ^Y in Eq. (45a), Eroshkin & Pashkevich (2007) in fact obtained ψ_so(t) instead of ${\psi }_{\text{SO}}^{*}\left(t\right)$.

In subsequent studies (e.g., Pashkevich 2016), the geodetic variations in Ecliptic Euler angles were described and took into account σ^X and σ^Y. Here, we compute the precession rate in ecliptic Euler angles for Mars, by multiplying 6.754 mas yr⁻¹ by $\left(\cos i_0-\sin i_0\cot {\epsilon }_0^{*}\cos \left({\psi }_0^{*}-{\text{Ω}}_0\right)\right)$, $-\sin i_0\sin \left({\psi }_0^{*}-{\text{Ω}}_0\right)$, and $\left(\sin i_0\cos \left({\psi }_0^{*}-{\text{Ω}}_0\right)\right)/\sin {\epsilon }_0^{*}$, respectively. With i₀ = 1°.84973, Ω₀ = 49°.5581, ${\psi }_0^{*}=82°.9071$, ${\epsilon }_0^{*}=26°.7179$, we find 6.389, −0.120 and 0.405 mas yr⁻¹ in ψ*‚ ε* and ϕ*, respectively. For the rate in obliquity and rotation angle, we obtain values that are consistent with Pashkevich (2016), but with opposite signs. The geodetic rate in longitude of Pashkevich (2016) is 7.114 mas yr⁻¹. We believe this value was obtained erroneously as a result of confusion with regard to the sign for the ecliptic obliquity, ${\epsilon }_0^{*}$.

Fig. 6 Signature in the Doppler observable of the relativistic corrections in the rotation angle ϕ: (a) Signature (in mm s−1) of the periodic relativistic variations $\left({\left[\varphi \right]}_{\text{GR}}\left(t\right)-{\dot{\varphi }}_{\text{GR}}t\right)$ in the rotation angle as a function of time, using the RISE timing (November 2018–April 2022). The blue envelope uses a simulated continuous timing. The gray dashed line represents the periodic relativistic variations in the rotation angle δϕ, arbitrarily rescaled. (b) Signature (in mm s−1) of both the linear and periodic relativistic terms [ϕ]GR(t) in the rotation angle, using the RISE timing.

Fig. 7 Signature in the Doppler observable of the relativistic corrections in the longitude angle ψ: (a) Signature in the Doppler observable (in μm s−1) of the small geodetic nutations in the longitude angle as a function of time, using the RISE timing (Nov. 2018–April 2022). The blue envelope uses a simulated continuous timing. The gray dashed line represents the geodetic nutations, arbitrarily rescaled. (b) Signature in the Doppler observable (in mm s−1) of all the geodetic terms (nutations and the linear term), using the RISE timing.

5.2 Rotation variations due to time coordinate transformation

In Table 6, we present estimates for the Sun contribution to the relativistic variations in ϕ, for a rotation model expressed in the BCRS, assuming that the planets follow Keplerian orbits (see toy model of Sect. 3.1.1). This two-body contribution tends to increase with increasing distance from the Sun, the maximum being reached for Saturn. The linear term includes for all planets the L_B contribution for the rescaling between TCB and TDB (see Eq. (28)).

Table 7 provides estimates for the contribution related to the Sun motion relative to the SSB, based on the toy model of Sect. 3.1.2. For each planet, seven terms at synodic periods corresponding the indirect effects of the other planets are computed. Most of these contributions have negligible amplitude. Only the giant planets induce indirect effects on the other planets larger than 0.01 mas in amplitude, and this effect is almost zero on Mercury and Venus. The largest term (with an amplitudes above 4 mas) applies to Saturn and is caused by Jupiter.

In Table 8, we compute the contribution related to the direct effect of the planets on each other, based on the toy model of Sect. 3.1.3. For each planet, we give the total rate of the linear term due to all the other planets and the amplitudes of the seven individual terms at synodic periods. The linear terms are very small, as already noticed for Mars. As for the indirect terms, the synodic direct terms are mainly induced by the giant planets, with quasi zero effect for Mercury and Venus, and the largest term (with an amplitudes above 1 mas) applies to Saturn and is caused by Jupiter.

The results presented in this section depend on the chosen values for the parameters and in particular on the values of the eccentricities. Here, we use the eccentricities obtained from the mean orbital elements of Simon et al. (2013). If necessary, the results can be refined by a fit to a numerical integration (as in Sect. 3.2) or from a semi-analytical solution (as in Sect. 3.3).

6 Discussion and conclusions

We have included relativistic corrections in the orientation and rotation model of Mars expressed in the BCRS. We estimated the corrections in the Euler angles (ψ, ε, ϕ) describing the orientation of a frame attached to the surface of Mars with respect to its mean orbit. Given the current accuracy on radioscience orbiter and lander data, a precision of 0.1 mas in the relativistic corrections is required to avoid errors in the interpretation of measurements of Mars rotation in terms of local physics. An accurate estimation of the relativistic corrections in rotation is also useful to define IAU standards for the rotation and orientation of Mars (Yseboodt et al., in prep.).

We first considered the relativistic terms that impact directly the rotation, finding that only the geodetic precession induces a significant effect and this is only in longitude, ψ. Then we investigated the terms that arise in the rotation angle, ϕ, because of the time coordinate transformation between a local Mars reference fame and the BCRS. For the longitude correction, our results are in agreement with previous findings, whereas this is not the case for the rotation correction. There is no significant relativistic correction that applies to the obliquity, ε.

Therefore, our recommendations for the relativistic corrections in Mars’ Euler angles (in mas) are as follows: $${\left[\psi \right]}_{\text{GR}}\left(t\right)=6.754t+0.565\sin l^{\prime },$$(46a) $${\left[\epsilon \right]}_{\text{GR}}\left(t\right)=0,$$(46b) $$\begin{array}{ll}{\left[\varphi \right]}_{\text{GR}}\left(t\right)\hfill & =7.3088d-166.954\sin l^{\prime }-7.783\sin 2l^{\prime }\hfill \\ \hfill & -0.544\sin 3l^{\prime }+0.567\sin \left(\frac{2\pi }{2.235294}t+320°.997\right)\hfill \\ \hfill & +0.102\sin \left(\frac{2\pi }{2.009124}t+303°.752,\right)\hfill \end{array}$$(46c)

with t and d the time in years and days, respectively, and l′ as the mean anomaly of Mars as given in Eq. (41c). For the longitude angle, we kept the linear and annual terms of Eq. (13), estimated from a Keplerian toy model, which is consistent with the results of Baland et al. (2020). The precision on those terms is of about 0.05%, as estimated from the difference with a semi-analytical derivation based on VSOP87 ephemerides (not shown here). The linear term in [ψ]_GR of 6.754 mas yr⁻¹ is important, making a difference in the angle up to 135 mas around the year 2020. For the rotation angle, we kept the linear term and the five largest periodic terms, based on an average of the fit to the numerical solution of Table 2 and of the semi-analytical solution of Table 3. The linear term in [ϕ]_GR of 7.3 mas day⁻¹ is very large, shifting the rotation angle by 53,000 mas around the year 2020, namely, moving the prime meridian by almost 850 m in 20 yr. The accuracy of the linear term is of the order of 10⁻⁵%, the difference between the two solutions. The precision on the periodic terms is better than 0.01 mas. The last two terms are at the Mars-Jupiter (2.24yr) and Mars-Saturn (2.01 yr) synodic periods. We note that the phase of these synodic terms is not obtained as the difference between the phases of the mean anomalies of Mars and of the other planet.

Our recommendation for the expression of [ϕ]_GR(t) replaces the estimate of Yoder & Standish (1997, reminded in Eq. (3)), where only the three main periodic terms are given, with an error of about 9 mas on the annual term. Such a difference can already have an effect in the radioscience data analysis, since the periods of these terms are the same as the periods of the rotation variations induced by atmosphere and surface dynamics. The synodic terms are here computed for the first time. Since their period is close the orbital period of Mars (1.88 yr), we recommend to include them in the a priori rotation model of Mars in order to avoid any contamination of the rotation amplitudes fitted to the radio-science data. Not taking them into account would likely affect the estimate of the annual term in ϕ by 0.6 mas and 0.1 mas, respectively (if the annual term absorbs their full contribution).

The applied methods (analytical, numerical, and semi-analytical) presented in this study can be extended to other bodies orbiting the Sun. In Sect. 5.2, we already demonstrated an application of the analytical toy model to other planets of the Solar System. These results can be refined by a fit to a numerical integration or from a semi-analytical solution. The methods should also be upgraded to moons, such as the Galilean moons.

Acknowledgements

We thank the editor and an anonymous reviewer for their careful reading and comments that have helped to improve our paper. This work was financially supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. This is InSight contribution ICN 304.

References

All Tables

All Figures

Fig. 1 Euler angles between the rotating body frame of Mars (axes XYZ) and the inertial frame associated with the mean orbit of Mars of 1980 (axes xyz, the x-axis is in the direction of the ascending node of Mars's orbit over the Earth ecliptic of epoch J2000). The X-axis of the BF is chosen as the prime meridian defined in the IAU convention (Archinal et al. 2018). The spin axis longitude ψ is measured from the x-axis to the autumn equinox, ϕ is measured from the equinox to the X-axis, and the obliquity, ε, is the angle from the z-axis to the Z-axis – or the inclination of the BF equator over the IF xy plane.

In the text

	Fig. 2 Difference between the [ϕ]GR(t) obtained by numerically integrating Eq. (17), using the DE440 planetary ephemerides and using Eq. (25) and the fitted series from Eq. (37). This curve provides an estimate of the accuracy of the fitted analytical model provided by the coefficients from Table 2.
In the text

	Fig. 3 Impact of the 1/c4 term from Eq. (17) on the BCRS modeling of Mars’s rotation (top). Impact of the Sun quadrupole moment $\left(J_2^{\odot }\right)$ on the BCRS modeling of Mars’s rotation (bottom).
In the text

	Fig. 4 Difference between the solutions for [ϕ]GR(t) obtained using VSOP87 and DE440 ephemerides via numerical integrations (left). Difference between the numerical integration and the semi-analytical solution using VSOP87, with only the 11 terms of Table 3 in blue, with all terms in black (right).
In the text

	Fig. 5 Difference between the [ϕ]GR(t) of the semi-analytical solution using VSOP87 (with only the 11 terms of Table 3) and the fit of the numerical integration with DE440 ephemerides.
In the text

Fig. 6 Signature in the Doppler observable of the relativistic corrections in the rotation angle ϕ: (a) Signature (in mm s−1) of the periodic relativistic variations $\left({\left[\varphi \right]}_{\text{GR}}\left(t\right)-{\dot{\varphi }}_{\text{GR}}t\right)$ in the rotation angle as a function of time, using the RISE timing (November 2018–April 2022). The blue envelope uses a simulated continuous timing. The gray dashed line represents the periodic relativistic variations in the rotation angle δϕ, arbitrarily rescaled. (b) Signature (in mm s−1) of both the linear and periodic relativistic terms [ϕ]GR(t) in the rotation angle, using the RISE timing.

In the text

Fig. 7 Signature in the Doppler observable of the relativistic corrections in the longitude angle ψ: (a) Signature in the Doppler observable (in μm s−1) of the small geodetic nutations in the longitude angle as a function of time, using the RISE timing (Nov. 2018–April 2022). The blue envelope uses a simulated continuous timing. The gray dashed line represents the geodetic nutations, arbitrarily rescaled. (b) Signature in the Doppler observable (in mm s−1) of all the geodetic terms (nutations and the linear term), using the RISE timing.

In the text