© The Authors 2025

1 Introduction

The existence of a molten core in Mercury has been inferred since the detection of its magnetic field by the Mariner 10 spacecraft (Ness et al. 1974). Shortly thereafter, Peale (1976) discussed methods to independently estimate the polar moments of inertia, C for the entire planet and C_m for the mantle alone, using the planet’s obliquity and forced libration, knowing the gravitational coefficients C₂₀ and C₂₂. This estimation is feasible due to the core-mantle coupling, which is weak enough for the core to avoid following the mantle’s rapid 88-day libration period, yet strong enough for the core to follow the mantle’s slow 327000-year nodal precession (Peale 1976)¹.

Subsequent measurements have been conducted using Earth-based radar observations for libration (Margot et al. 2007) and data from the MESSENGER spacecraft for obliquity and gravity coefficients (Mazarico et al. 2014). These observations nicely confirm the predictions of Peale (1976), with current estimates placing C/(MR2) at 0.333 ± 0.005 and C_m/C at 0.443 ± 0.019 (Genova et al. 2019; see also Goossens et al. 2022).

With this value of the mantle’s moment of inertia, the period of Mercury’s free oscillation in longitude is close to Jupiter’s 12-year orbital period (Margot et al. 2007). The amplitude and phase of the forced libration at this period should thus allow for the refinement of the mantle’s moment of inertia (Peale et al. 2009). However, observations over a sufficiently long timespan are still missing to properly fit these two quantities (Margot et al. 2007).

Other physical properties that can be probed by rotation include tidal dissipation and core-mantle interactions, such as viscous coupling, topographic coupling, electromagnetic coupling, and pressure coupling (Peale et al. 2014). The pressure torque, related to the flattening of the core, is a conservative torque that mainly modifies the equilibrium obliquity (Peale et al. 2014).

All other dissipative effects influence the amplitude and phase of the forced libration at 12 years (Peale et al. 2009), and they also displace the spin axis from the plane containing the normal to the invariant plane and the orbit pole, known as the Cassini plane (Peale et al. 2014).

The dissipative effects on Mercury’s spin axis orientation have been analysed numerically by Peale et al. (2014), and tides in particular, characterised by the ratio k₂/Q of the Love number k₂ divided by the quality factor Q, have been studied analytically for a solid Mercury without a molten core by Baland et al. (2017). Notably, all core-mantle dissipative mechanisms have quite similar signatures (Peale et al. 2014) and a significantly stronger effect than tides (Peale 2005). The addition of a solid inner core has been studied both numerically and analytically by (Peale et al. 2016; Dumberry 2021; MacPherson & Dumberry 2022).

Concurrently, observational studies of Earth-based measurements and MESSENGER data have constrained Mercury’s spin axis orientation, suggesting that Mercury’s spin indeed deviates from the Cassini plane by a few arcseconds (Margot et al. 2012; Mazarico et al. 2014; Stark et al. 2015a; Stark et al. 2015b; Verma & Margot 2016; Genova et al. 2019; Konopliv et al. 2020; Bertone et al. 2021). However, current observational uncertainties still leave room for further refinement, particularly in constraining the core-mantle interactions.

The upcoming BepiColombo mission, a joint endeavour between the European Space Agency (ESA) and the Japan Aerospace Exploration Agency (JAXA), launched on October 19, 2018, represents a significant advancement in our ability to probe Mercury’s geophysical properties. Upon its scheduled orbital insertion in 2026, BepiColombo will deploy its suite of state-of-the-art instruments to achieve unprecedented accuracy in measuring geodetic parameters such as the gravitational coefficients C₂₀ and C₂₂, the amplitude of longitudinal librations, the orientation of the spin axis, and tidal deformation (Genova et al. 2021). These improved measurements will provide a more precise assessment of Mercury’s rotational motion, offering constraints on its interior structure (Steinbrügge et al. 2021).

In this study, we revisit Mercury’s rotational dynamics and derive analytical formulae for the orientation of its spin axis, extending the model presented in Baland et al. (2017) to include a fluid core. These formulae are designed to avoid the need for numerical integrations as in Peale et al. (2014). Given that the different core-mantle dissipative processes discussed by Peale et al. (2014) have similar signatures, here, we only consider the viscous coupling at the core-mantle boundary (CMB), which is proportional to the relative angular velocity between the core and the mantle (e.g. Peale et al. 2007).

While acknowledging the comprehensive model of MacPherson & Dumberry (2022) that includes the solid core, our initial focus is on refining the formula for the deviation of the spin axis from the Cassini plane induced by tides, addressing an inconsistency identified in Baland et al. (2017) and subsequently reproduced in MacPherson & Dumberry (2022). Additionally, we highlight a non-linear effect that emerges at low viscosity for specific values of the pressure torque, associated with the substantial tilt of the core’s angular velocity as depicted in Peale et al. (2014, Fig. 3) and discussed in Boué (2020). The exploration of more realistic interior structures is reserved for future work.

Our analysis employs the Lagrangian formalism developed in Ragazzo & Ruiz (2015) and Boué (2020), which is mathematically equivalent to the angular momentum approach followed by Baland et al. (2017). In Sections 2 and 3, we pursue a similar series of approximations as in Baland et al. (2017) to obtain analytical solutions. Section 2 is dedicated to the rigid body approximation while the tidal deformation is accounted in Section 3. The key distinction, beyond the incorporation of the fluid core in Section 4, lies in our consideration of the complete 3D rotation matrix defining Mercury’s orientation (three degrees of freedom), whereas Baland et al. (2017) focused solely on the spin axis (two degrees of freedom). As demonstrated in our results in Section 5, this methodological difference significantly impacts the outcomes, despite both approaches yielding consistent results when considering only the projection of the torque orthogonal to the spin axis. Our conclusions are presented in Section 6.

In this paper, we follow the notation of Boué (2020) and represent each vector v by a skew-symmetric matrix v such that $$\widehat{\mathit{v}}=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\to \mathit{v}=\left(\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right).$$(1)

The dot product 〈, 〉 of two matrices is defined as $$\langle \mathit{a},\mathit{b}\rangle \equiv \frac{1}{2}Tr\left({\mathit{a}}^T\mathit{b}\right)=\widehat{a}\cdot \widehat{b},$$(2)

where Tr denotes the trace operator, and the norm $\sqrt{\langle \mathit{v},\mathit{v}\rangle }$ of a vector v is denoted ∥v∥.

2 Mercury as a rigid body

Here, we retrieve the results of Baland et al. (2017) using the formalism of Boué (2020). The goal is to introduce the methodology on a simple case before adding complexity, namely deformation and a liquid core.

For the sake of generality, hereafter we refer to Mercury simply as a body ℬ of mass M and radius R orbiting a central point mass object m₀. Let R_ℬ be the rotation matrix between an inertial reference frame ℱ_inv and a frame ℱ_ℬ attached to ℬ. For the inertial frame, we take a basis (i_inv, j_inv, k_inv) whose third axis is normal to the invariant plane with respect to which the orbit is assumed to have a constant inclination i. As for the body fixed frame, we choose the basis composed of the principal axes of inertia (I, J, K), in which the moments of inertia are denoted A, B and C along the first, second and third axis, respectively.

Let Ω be the longitude of the ascending node of the orbit reckoned from the invariant plane. We assume the planet to be in a p : 1 spin-orbit resonance, meaning that its angular velocity with respect to the inertial frame is $$\begin{array}{cc}{\mathit{\omega }}_{ℬ}=\dot{\text{Ω}}{\mathit{k}}_{\text{inv}}+{\omega }_p\mathit{K},& {\omega }_p\equiv pn+\dot{\varpi }-\end{array}\dot{\text{Ω}},$$(3)

with $n=\sqrt{G{m}_0/a^3}$ the orbital mean motion and ϖ the longitude of the periapsis. In practice, for Mercury, we have p = 3/2. In the expression of the mean motion, G is the gravitational constant and a the semi-major axis.

We consider that the position of the Cassini state corresponds to a fixed orientation of the vector K in a frame rotating at the orbital nodal precession frequency. Note that in the literature, Cassini states are sometimes defined using the angular momentum vector (Noyelles et al. 2011) instead of the body polar axis. However, both axes are very close, and for our purpose the two definitions are equivalent. Let us then introduce the frame ℱ_orb attached to the orbit whose basis (i, j, k) is such that i is along the orbital ascending node with respect to the invariant plane and k is along the orbital angular momentum. The transformation of coordinates between the orbital frame ℱ_orb and the inertial invariant frame ℱ_inv is given by R_inv = R₃(Ω) R₁ (i) where R₁ denotes a rotation around the first axis and R₃ around the third axis². The rotation matrix between the body fixed frame ℱ_ℬ and the inertial frame ℱ_inv can itself be decomposed as follows $${\mathit{R}}_{ℬ}={\mathit{R}}_3\left(\text{Ω}\right){\mathit{R}}_1\left(i\right)\mathit{R}{\mathit{R}}_3\left(\theta p\right),$$(4)

with θ_p = pM + ϖ − Ω + γ₀ the rotation angle of the planet, γ₀ an offset being equal to 0 or π/2 depending on the resonance (Wisdom et al. 1984, page 139), M the mean anomaly, and R the rotation defining the orientation of the Cassini states. Were the body not in a Cassini state, the expression of R_ℬ would still hold, but the angular velocity would have an extra term $\mathit{\omega }=\dot{\mathit{R}}{\mathit{R}}^T$ accounting for the time derivative of R, namely $${\mathit{\omega }}_{ℬ}=\dot{\text{Ω}}{\mathit{k}}_{\text{inv}}+\mathit{\omega }+{\omega }_p\mathit{K}.$$(5)

In the following, we provide the equations of motion in terms of the variables (ω, R), the Cassini state being an equilibrium of the secular problem (averaged over the orbital mean anomaly) corresponding to ω = 0 and R = constant. When deformation is taken into account, the matrix of inertia (or equivalently the traceless symmetric B-matrix, as defined in Ragazzo & Ruiz (2015)) is allowed to vary, and an equation of motion is added for this variable.

The Lagrangian describing the motion of (ω, R), written in the orbital frame ℱ_orb is (e.g. Boué 2020) $$\mathit{ℒ}=\frac{{\text{I}}_{\circ }}{2}\left({\Vert {\mathit{\omega }}_{ℬ}\Vert }^2-2\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\mathit{\omega }}_{ℬ}^2\rangle \right)+3{\text{I}}_{\circ }\frac{G{m}_0}{a^3}\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,\mathit{S}\rangle ,$$(6)

with I_∘ ≡ (A + B + C)/3 the mean moment of inertia of the body ℬ, B′ = R₃(θp)BR₃(–θ_p), the B-matrix written in the preces-sional frame as defined in Ragazzo et al. (2022, Sect. 7) such that RB′R^T is expressed in the orbital frame ℱ_orb, and S a traceless symmetric matrix encapsulating the orbital dependency of the spin-orbit potential energy. We only consider Mercury’s gravity field at second order, as the higher orders measured by MESSENGER (Smith et al. 2012; Mazarico et al. 2014) have negligible effects on the orientation of the spin axis (Noyelles & Lhotka 2013; Baland et al. 2017). We decompose B and S onto real spherical harmonic matrices as defined in Boué (2020), namely, $$\mathit{B}=\frac{2}{3}\frac{M{R}^2}{r^3}\left(C_{20}{\mathit{Y}}_{20}+C_{22}{\mathit{Y}}_{22}\right),$$(7) $$\mathit{S}=\frac{a^3}{r^3}\left(-\frac{1}{2}{\mathit{Y}}_{20}+\frac{1}{6}\left(\cos \left(2\mathit{v}+2\varpi -2\text{Ω}\right){\mathit{Y}}_{22}+\sin \left(2\mathit{v}+2\varpi -2\text{Ω}\right){\mathit{Y}}_{2-2}\right)\right),$$(8)

where r is the orbital separation and v the true anomaly. Parameters C₂₀ and C₂₂ are Stokes coefficients related to the moments of inertia according to $$C_{20}\equiv \frac{A+B-2C}{2M{R}^2},$$(9) $$C_{22}\equiv \frac{B-A}{4M{R}^2}.$$(10)

The equations of motion are (Boué 2020) $$\begin{array}{cc}\frac{\text{d}}{\text{d}t}\frac{\partial \mathit{𝓛}}{\partial \mathit{\omega }}=\left[\mathit{\omega },\frac{\partial \mathit{𝓛}}{\partial \mathit{\omega }}\right]+\widehat{\mathit{𝓙}}\left(\mathit{𝓛}\right),& \dot{\mathit{R}}=\mathit{\omega }\mathit{R}\end{array},$$(11)

where $\widehat{𝓙}$ is the spin operator defined for all function f (R) as $\widehat{𝓙}\left(f\right)\equiv \left(\partial f/\partial \mathit{R}\right){\mathit{R}}^T$ and [, ] is the usual commutator of two matrices : [A, B] ≡ AB – BA. The equilibrium is reached when the body’s angular momentum $$\mathit{\pi }\equiv \frac{\partial 𝓛}{\partial \mathit{\omega }}={\text{I}}_{\circ }\left({\mathit{\omega }}_{ℬ}+\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T{\mathit{\omega }}_{ℬ}+{\mathit{\omega }}_{ℬ}\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T\right)$$(12)

averaged over an orbital period remains constant for ω = 0. Substituting the expression of the Lagrangian in the equation of motion, one gets $$\dot{\mathit{\pi }}=-\dot{\text{Ω}}\left[{\mathit{k}}_{\text{inv}},\mathit{\pi }\right]+3{\text{I}}_{\circ }{n}^2\left[\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,\mathit{S}\right].$$(13)

Setting the time average of the body’s angular momentum rate $\dot{\mathit{\pi }}$ equal to 0, we obtain an equation whose unknown is the rotation matrix R defining the orientation of the Cassini state in the orbital frame ℱ_orb. In the case of Mercury, the deviation of the figure-axis K from the orbit normal k is very small (about 2 arcmin). It is thus legitimate to linearise Eq. (13) over the tilt angle ε = ε₁ i + ε₂ j + ε₃ k, where ε₁ represents the in-plane offset of the Cassini state, ε₂ the deviation from the Cassini plane, and ε₃ the rotation lag of the body around its spin axis or, equivalently, the centre of libration in longitude. With this approximation, R ≈ 𝕀+ ε where 𝕀 is the identity matrix and RB′R^T ≈ B′ + [ε, B′]. Moreover, the precession rate of the orbit $\dot{\text{Ω}}$ is much smaller than the orbital mean motion n in absolute value ($\left|\dot{\text{Ω}}/n\right|\approx 7\times {10}^{-7}$), therefore we can safely ignore terms in ${\dot{\text{Ω}}}^2$ before those in n². Applying the two simplifications in the expression of $\dot{\mathit{\pi }}$ (Eq. (13)) and averaging over the mean anomaly, one gets $$\begin{array}{l}\left[\begin{array}{ccc}D& 0& 0\\ 0& D& 0\\ 0& N& -4{\kappa }_{22}\end{array}\right]\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\end{array}\right)=-\left(\begin{array}{c}N\\ 0\\ 0\end{array}\right)\\ +{\kappa }_{\mathit{\omega }}\left(\begin{array}{ccc}\cos \left(2\varpi -2\text{Ω}\right)& \sin \left(2\varpi -2\text{Ω}\right)& 0\\ \sin \left(2\varpi -2\text{Ω}\right)& -\cos \left(2\varpi -2\text{Ω}\right)& 0\\ 0& 0& 0\end{array}\right)\end{array}$$(14)

with $$D=\kappa +C{\omega }_p\dot{\text{Ω}}\cos i,$$(15)

and $$N=\dot{\text{Ω}}\sin i\left(C{\omega }_p-M{R}^2{C}_{20}\dot{\text{Ω}}\cos i\right).$$(16)

The polar moment of inertia is given by $C={\text{I}}_{\circ }-\frac{2}{3}M{R}^2{C}_{20}$. In the above expressions, we applied a notation similar to that introduced in Baland et al. (2017) up to a factor $\frac{3}{2}n$. More explicitly, we have $${\kappa }_{20}=-\frac{3}{2}{n}^{2}M{R}^2{C}_{20}{X}_0^{-3,0}\left(e\right),$$(17) $${\kappa }_{22}=3\sigma {}_p{n}^{2}M{R}^2{C}_{22}{X}_{2p}^{-3,2}\left(e\right),$$(18)

and $${\kappa }_{\mathit{\omega }}=3\sigma {}_p{n}^{2}M{R}^2{C}_{22}{X}_{2p}^{-3,0}\left(e\right),$$(19)

where the $X_k^{n,m}\left(e\right)$ are Hansen coefficients of the eccentricity e defined such that $${\left(\frac{r}{a}\right)}^n\exp \left(imv\right)\equiv {\displaystyle \sum _{\mathit{k}=-\infty }^{\infty }{X}_k^{n.m}}\left(e\right)\exp \left(ikM\right),$$(20)

and ${\sigma }_p\equiv \cos \left(2{\gamma }_0\right)=sign\left(X_{2p}^{-3,2}\left(e\right)\right)$ (Wisdom et al. 1984). The system of equations (14) was obtained with the help of the computer algebra system TRIP (Gastineau & Laskar 2011). This system includes trigonometric terms in 2ϖ – 2Ω due to the apsidal precession of the orbit as in Baland et al. (2017, Eqs. (34, 35)). These terms are responsible for a nutation of the spin axis at a slow period of about 84 kyr (Baland et al. 2017). To solve Eq. (14), we proceeded as in Baland et al. (2017), i.e. we first neglected the contribution due to the apsidal precession of the orbit, which amounts to set κ_ω ~ e^2p equal to zero. In that case, we get a first order solution ${\epsilon }_1^{\left(1\right)}={\epsilon }_{\text{Ω}}$ and ${\epsilon }_2^{\left(1\right)}={\epsilon }_2^{\left(1\right)}=0$ with $${\epsilon }_{\text{Ω}}=-\frac{N}{D}=-\frac{\dot{\text{Ω}}\sin i\left(C{\omega }_p-M{R}^2{C}_{20}\dot{\text{Ω}}\cos i\right)}{\kappa +C{\omega }_p\dot{\text{Ω}}\cos i}.$$(21)

In a second step, we inserted this first order solution in the right-hand side of Eq. (14) and solved again the system. The result is $${\epsilon }_1={\epsilon }_{\omega }\cos \left(2\varpi -2\text{Ω}\right)+{\epsilon }_{\text{Ω}},$$(22a) $${\epsilon }_2={\epsilon }_{\omega }\sin \left(2\varpi -2\text{Ω}\right)+{\epsilon }_{\text{Ω}},$$(22b) $${\epsilon }_3=\frac{N}{4{\kappa }_{22}}{\epsilon }_{\omega }\sin \left(2\varpi -2\text{Ω}\right)+{\epsilon }_{\text{Ω}},$$(22c)

with $${\epsilon }_{\omega }=\frac{{\epsilon }_{\text{Ω}}{\kappa }_{\mathit{\omega }}}{\kappa +C{\mathit{\omega }}_p\dot{\text{Ω}}\cos i}.$$(23)

The results of this calculation from Eq. (21) to Eq. (23) are very similar to those in Baland et al. (2017). The main difference is the libration angle ε₃ that is absent in Baland et al. (2017).

Let us rewrite the relevant formulae neglecting the nodal precession frequency $\left|\dot{\text{Ω}}\right|$ before the mean motion n in a simpler form : $${ϵ}_{\text{Ω}}\approx \frac{C}{M{R}^2}\frac{p\left|\dot{\text{Ω}}\right|}{n}\frac{\sin i}{-\frac{3}{2}{C}_{20}{X}_0^{-3,0}\left(e\right)+3{\sigma }_p{C}_{22}{X}_{2p}^{-3,2}\left(e\right)},$$(24) $${ϵ}_{\omega }\approx \frac{2{\sigma }_p{C}_{22}{X}_{2p}^{-3,0}\left(e\right)}{-C_{20}{X}_0^{-3,0}\left(e\right)+2{\sigma }_p{C}_{22}{X}_{2p}^{-3,2}\left(e\right)}{\epsilon }_{\text{Ω}}.$$(25)

3 Mercury as a deformable body

Here, we consider the case where the body is deformable. As in Ragazzo & Ruiz (2015), we introduce an effective macroscopic rigidity μ accounting for the elastic response of the body as a whole according to the elastic potential energy $$U_{\text{el}}=\frac{{\text{I}}_{\circ }}{2}\mu {\Vert \mathit{B}-{\mathit{B}}_0\Vert }^2,$$(26)

where B₀ is a permanent deformation induced by a prestress (Ragazzo et al. 2022). The effective rigidity μ is related to the Love number k₂ and to the tidal parameter (Baland et al. 2017, Eq. (B.12)) $$q_t=-3\frac{G{m}_0}{GM}{\left(\frac{R}{a}\right)}^3$$(27)

through $$\frac{1}{\mu }=-\frac{k_2{q}_t}{9n^2}\frac{M{R}^2}{{\text{I}}_{\circ }}.$$(28)

For a smooth comparison with Baland et al. (2017), we also introduce the ratio of the centrifugal acceleration to the gravitational acceleration $q_r=-\frac{1}{3}{p}^2{q}_t$. With the elastic potential energy, the Lagrangian of the problem now reads $$\begin{array}{l}𝓛=\frac{{\text{I}}_{\circ }}{2}\left({\Vert {\mathit{\omega }}_{ℬ}\Vert }^2+2\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\mathit{\omega }}_{ℬ}{\mathit{\omega }}_{ℬ}^T\rangle \right)-\frac{{\text{I}}_{\circ }}{2}\mu {\Vert \mathit{B}-{\mathit{B}}_0\Vert }^2\\ +3{\text{I}}_{\circ }\frac{G{m}_0}{a^3}\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,\mathit{S}\rangle .\end{array}$$(29)

The viscous deformation of the body as a whole is controlled by a Rayleigh dissipation function 𝒟 that is a function of an effective macroscopic viscosity η such that $$\mathit{𝒟}=\frac{{\text{I}}_{\circ }}{2}\eta {\Vert \dot{\mathit{B}}\Vert }^2.$$(30)

From the effective macroscopic viscosity η and the effective macroscopic rigidity μ, one can define a Maxwell relaxation time τ_M = η/μ that is related to the tidal quality factor Q = 1/(nτ_M). Introducing the phase shift ζ ≡ 1/Q = nτ_M as in Baland et al. (2017), we get the relation $$\eta =-\frac{9n\zeta }{k_2{q}_t}\frac{{\text{I}}_{\circ }}{M{R}^2}.$$(31)

The equations of motion (11) for (ω, R) still hold. We now add the equation describing the evolution of the B-matrix, namely, $$\frac{\text{d}}{\text{d}t}\frac{\partial 𝓛}{\partial \dot{\mathit{B}}}=\frac{\partial 𝓛}{\partial \mathit{B}}-\frac{\partial \mathit{𝒟}}{\partial \dot{\mathit{B}}}.$$(32)

Since we discard the kinetic energy associated with the deformation of the body (there is no Ḃ dependency in the Lagrangian), the equation of motion (32) is equal to 0. This approximation is justified in Correia et al. (2018). Let δB ≡ B – B₀ be the deviation of the B-matrix from its relaxed state B₀. From the expression of the Lagrangian (29) and of Rayleigh dissipation function (30), the equation of motion (32) reads $$\delta \mathit{B}+\frac{\eta }{\mu }\delta \dot{\mathit{B}}=\frac{1}{\mu }\mathit{F}$$(33)

with F = F_r + F_t and $${\mathit{F}}_r={\mathit{R}}_3\left(-{\theta }_p\right){\mathit{R}}^T\left({\mathit{\omega }}_{ℬ}{\mathit{\omega }}_{ℬ}^T\right)\mathit{R}{\mathit{R}}_3\left({\theta }_p\right)-\frac{2}{3}{\Vert {\mathit{\omega }}_{ℬ}\Vert }^2\text{I},$$(34) $${\mathit{F}}_t=3n^2{\mathit{R}}_3\left(-{\theta }_p\right){\mathit{R}}^T\mathit{S}\mathit{R}{\mathit{R}}_3\left({\theta }_p\right).$$(35)

Given that the influence of deformation on the spin axis orientation is weak, only the leading-order terms are retained in the evaluation of the ‘force’ F on the left-hand side of Eq. (33). Consequently, we disregard the obliquity, assuming the spin axis to be perpendicular to the orbital plane, and neglect the precession frequency relative to the spin rate. With this approximation, $${\mathit{F}}_r\approx {\mathit{\omega }}_p^2\left(\mathit{k}{\mathit{k}}^T-\frac{2}{3}\text{𝕀}\right),$$(36) $${\mathit{F}}_t\approx 3n^2{\mathit{R}}_3\left(-{\theta }_p\right)\mathit{S}{\mathit{R}}_3\left({\theta }_p\right).$$(37)

The static part $F^{\text{stat}}\equiv \overline{F}$ of F generates a constant deformation δB^stat for which the Love number k₂ admits the fluid limit k_f. We get $$\delta {\mathit{B}}^{\text{stat}}=\frac{2}{3}\left(\delta C_{20}^{\text{stat}}{\mathit{Y}}_{20}+\delta C_{22}^{\text{stat}}{\mathit{Y}}_{22}\right)$$(38)

with $$\delta C_{20}^{\text{stat}}=\frac{k_f}{6}\left(-2q_r+q_t{X}_0^{-3,0}\left(e\right)\right),$$(39) $$\delta C_{22}^{\text{stat}}=\frac{k_f}{12}{q}_t{X}_0^{-3,0}\left(e\right)\cos \left(2{\gamma }_0\right).$$(40)

These two formulae are in perfect agreement with Baland et al. (2017, Eqs. (B.13–B.14)). Putting numbers into these expressions show that the hydrostatic parts $\delta C_{20}^{\text{stat}}$ and $\delta C_{22}^{\text{stat}}$ represent about 1% of the observed gravity field coefficients C₂₀ and C₂₂, respectively.

We now turn to the oscillatory part $\mathit{F}-\overline{\mathit{F}}$ that produces a significant out of plane offset of the Cassini state of about 1 arcsec according to Baland et al. (2017). At first order in τ_Μ = η/μ, the oscillatory (or periodic) part δB^peri of δB is $$\delta {\mathit{B}}^{\text{peri}}=\frac{1}{\mu }\left(\mathit{F}-\overline{\mathit{F}}-{\tau }_m\dot{\mathit{F}}\right),$$(41)

with μ evaluated at k₂ ≠ k_f. In our approximation, the centrifugal force remains static in the body frame. Consequently, we can compute the periodic component of δB using Eq. (41), substituting the force F with F_t.

Once δB^peri is known, the system (14) is recomputed with B′ replaced by B′ + R₃(θ_p)δB^periR₃(–θ_p) in Eqs. (12 and 13). The result is averaged over the mean anomaly and over the longitude of the periapsis to focus on the out of plane offset of the Cassini state. The calculation done with TRIP gives $$\left[\begin{array}{ccc}𝒟& 0& 0\\ 0& 𝒟& 0\\ 0& N& -4{\kappa }_{22}\end{array}\right]\left(\begin{array}{c}\delta {\epsilon }_1\\ \delta {\epsilon }_2\\ \delta {\epsilon }_3\end{array}\right)=-\left(\begin{array}{c}{\kappa }_{\zeta ,c}{\epsilon }_{\text{Ω}}\\ {\kappa }_{\zeta ,s}{\epsilon }_{\text{Ω}}\\ {\kappa }_{\zeta ,s}{\kappa }_{\zeta ,n}\end{array}\right),$$(42)

where δε = (δε₁, δε₂, δε₃) is a correction to the solution ε = (ε₁, ε₂, ε₃) obtained in the rigid case (Eq. (22)). The constants in this new system are defined as $${\kappa }_{\zeta ,c}\equiv \frac{1}{4}{k}_2{q}_t{n}^{2}M{R}^2\left({\left(X_{2p}^{-3,2}\left(e\right)\right)}^2+{\left(X_0^{-3,0}\left(e\right)\right)}^2\right),$$(43) $${\kappa }_{\zeta ,s}\equiv \frac{1}{2}{k}_2{q}_t\zeta n{\omega }_{p}M{R}^2{X}_0^{-6,0}\left(e\right),$$(44) $${\kappa }_{\zeta ,n}\equiv \frac{1}{2}{k}_2{q}_t\zeta nM{R}^2\left({\omega }_p{X}_0^{-6,0}\left(e\right)-2n\sqrt{1-e^2}{X}_0^{-8,0}\left(e\right)\right).$$(45)

The first line of Eq. (42) provides a small correction to the in-plane offset of the Cassini state ε_Ω of the order of 0.1% that becomes $${\epsilon }_{\text{Ω}}^{k_2}={\epsilon }_{\text{Ω}}\left(1-\frac{{\kappa }_{\zeta ,c}}{𝒟}\right).$$(46)

The third line corresponds to a rotation of the body around its figure axis of about 3 arcsec, and the second line gives the out of plane offset of the Cassini state, namely, δε₂ = ε_ζ with $${\epsilon }_{\zeta }=-\frac{{\kappa }_{\zeta ,s}{\epsilon }_{\text{Ω}}}{D}=-\frac{{\kappa }_{\zeta ,s}{\epsilon }_{\text{Ω}}}{\kappa +C{\omega }_p\dot{\text{Ω}}\cos i}.$$(47)

This result shows some similarities with Baland et al. (2017, Eq. (66)). Our denominator is directly expressed in terms of the entries of the B-matrix while in Baland et al. (2017) the same denominator (up to a coefficient cos i) explicitly shows the contribution of the static hydrostatic adjustment (related to the δB^stat-matrix) with respect to a situation where only the relaxed gravity field coefficients (related to the B₀-matrix) would have been taken into account. But since observations cannot distinguish these two constant components, we chose not to split the entries of the B-matrix (which means that in our notation, C, C₂₀ and C₂₂ do include the permanent hydrostatic deformation).

The primary difference between our result (47) and Baland et al. (2017, Eq. (66)) lies in the numerator. Our expression includes only the first term of Baland et al.’s numerator, omitting the cos i factor. To derive their Eq. (66), Baland et al. (2017) solved two components of their 3D angular momentum equation (Eq. (59)), where the unknown is the spin axis (ŝ(t) in their notation). The challenge arises because their Eq. (59) for dŝ/dt is not orthogonal to ŝ, despite ŝ being a unit vector. By definition, if ŝ ∙ ŝ = 1, then $\frac{\text{d}\widehat{s}}{\text{d}t}\cdot \widehat{s}$ should be zero.

To address this issue, Baland et al. (2017) would have considered the two components of the torque orthogonal to ŝ (or to the orbit pole k, given the very small obliquity of the Cassini state ε_Ω). They instead chose the two components in the invariant plane, resulting in their numerator containing the second component of the right-hand side of our Eq. (42) multiplied by cos i, plus the third component of the right-hand side of our Eq. (42) multiplied by sin i. As a result, the reasoning used in Baland et al. (2017) to obtain their final expression Eq. (66) is not correct. We nevertheless verified that our expression for ε_ζ can be derived from Baland et al.’s equation of motion (Eq. (59)) using the correct projections.

4 Influence of a fluid core

In this section, we take into account the presence of a molten core in the interior of Mercury. The mantle surrounding the core is still subject to tidal deformation. The fluid core is parameterised by its mean moment of inertia I_∘c and by the B-matrix $${\mathit{B}}_{\text{c}}\equiv \frac{2}{3}\frac{M{R}^2}{{\text{I}}_{\circ \text{c}}}\left(C_{20}^{\text{c}}{\mathit{Y}}_{20}+C_{22}^{\text{c}}{\mathit{Y}}_{22}\right),$$(48)

where the Stokes coefficients $C_{20}^{\text{c}}$ and $C_{22}^{\text{c}}$ are normalised by the product MR² of the total mass M and the full radius squared R² of the body. The matrix B_c implicitly defines the shape of the core mantle boundary (CMB). The mantle is equiva-lently characterised by its mean moment of inertia I_∘m and the B-matrix $${\mathit{B}}_{\text{m}}\equiv \frac{2}{3}\frac{M{R}^2}{{\text{I}}_{\circ m}}\left(C_{20}^{\text{m}}{\mathit{Y}}_{20}+C_{22}^{\text{m}}{\mathit{Y}}_{22}\right),$$(49)

where $C_{20}^{\text{m}}$ and $C_{22}^{\text{m}}$ also are normalised by MR². The overall body has a mean moment of inertia I_∘ = I_∘c + I_∘m and its gravity field is represented by the matrix B such that I_∘B = I_∘cB_c + I_∘mB_m. These three matrices rotate in concert with the rotation matrix of the mantle R. Under tidal perturbation, B_c and B_m could vary independently. We should expect two normal modes of vibration, a symmetric one where the two matrices oscillate in phase and another antisymmetric where they oscillate in antiphase. For the sake of simplicity, we only consider the symmetric mode and assume that each of these B-matrices remain proportional to the total B-matrix with a constant coefficient that we denote f_c for the core and f_m for the mantle, namely, $$\begin{array}{ccc}{\mathit{B}}_{\text{c}}\equiv f_c\mathit{B},& and& {\mathit{B}}_{\text{m}}\equiv f_m\mathit{B}.\end{array}$$(50)

By definition of the matrix B, the two coefficients f_c and f_m are related to each other by the formula $${\text{I}}_{\circ c}{f}_c+{\text{I}}_{\circ m}{f}_m={\text{I}}_{\circ }.$$(51)

The rotational state of the core is given by its angular velocity $${\mathit{\omega }}_{\mathit{B}}^{\text{c}}=\dot{\text{Ω}}{\mathit{k}}_{\text{inv}}+{\mathit{\omega }}_{\text{c}}$$(52)

where ω_c is the rotation speed with respect to the orbital frame. For the mantle, we chose the same decomposition as in the single layer model above, i.e. we define its angular velocity as $${\mathit{\omega }}_{\mathit{B}}^{\text{m}}=\dot{\text{Ω}}{\mathit{k}}_{\text{inv}}+{\mathit{\omega }}_p\mathit{K}.$$(53)

The variables of the problem are (ω_c, ω_m, R, B). At a Cassini state, the system is in a relative equilibrium characterised by ω_m = 0, ω_c = constant and R = constant.

The Lagrangian of the system is the kinetic energy of the core plus the kinetic energy of the mantle minus the spinorbit potential energy minus the elastic potential energy of deformation $$\begin{array}{l}𝓛=\frac{{\text{I}}_{\circ c}}{2}\left({\Vert {\mathit{\omega }}_{ℬ}^{\text{m}}\Vert }^2+2f_c\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\mathit{\omega }}_{ℬ}^{\text{c}}{\mathit{\omega }}_{ℬ}^{cT}\rangle \right)\\ +\frac{{\text{I}}_{\circ c}}{2}\left({\Vert {\mathit{\omega }}_{ℬ}^{\text{m}}\Vert }^2+2f_m\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\mathit{\omega }}_{ℬ}^{\text{m}}{\mathit{\omega }}_{ℬ}^{mT}\rangle \right)\\ +3{\text{I}}_{\circ }{n}^2\langle \mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,\mathit{S}\rangle -\frac{{\text{I}}_{\circ }}{2}\mu {\Vert \mathit{B}-{\mathit{B}}_0\Vert }^2.\end{array}$$(54)

To this Lagrangian, we added a Rayleigh dissipation function similar to the one used above in the single layer approach but with an extra term accounting for a viscous friction at the CMB, namely, $$𝒟=\frac{{\text{I}}_{\circ }}{2}\eta {\Vert \dot{\mathit{B}}\Vert }^2+\frac{{\text{I}}_{\circ }}{2}\nu {\Vert {\mathit{\omega }}_{ℬ}^{\text{c}}-{\mathit{\omega }}_{ℬ}^{\text{m}}\Vert }^2,$$(55)

where ν (not to be confused with the true anomaly v) has the dimension of a frequency related to the inverse of the characteristic timescale $${\tau }_{\nu }=\left(1-\frac{C_m}{C}\right)\frac{C_m}{C}\frac{1}{\nu }$$

required to damp the relative angular velocity of the core with respect to that of the mantle (Peale et al. 2009). In this formula, $C_m={\text{I}}_{\circ m}-\frac{2}{3}{f}_{m}M{R}^2{C}_{20}$ is the mantle’s polar moment of inertia. The equations of motion are $${\dot{\mathit{\pi }}}_{\text{c}}=\left[{\mathit{\omega }}_{\text{c}},{\mathit{\pi }}_{\text{c}}\right]-\frac{\partial 𝒟}{\partial {\mathit{\omega }}_{\text{c}}},$$(56a) $${\dot{\mathit{\pi }}}_{\text{m}}=\left[{\mathit{\omega }}_{\text{m}},{\mathit{\pi }}_{\text{m}}\right]-\frac{\partial 𝒟}{\partial {\mathit{\omega }}_{\text{m}}}-\mathit{𝓙}\left(\mathit{𝓛}\right),$$(56b) $$0=\frac{\partial \mathit{𝓛}}{\partial \mathit{B}}-\frac{\partial 𝒟}{\partial \dot{\mathit{B}}},$$(56c)

where the angular momenta are given by $${\mathit{\pi }}_{\text{c}}\equiv \frac{\partial 𝓛}{\partial {\mathit{\omega }}_{\text{c}}}={\text{I}}_{\circ ,c}\left({\mathit{\omega }}_{\mathit{ℬ}}^{\text{c}}+f_c\left(\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T{\mathit{\omega }}_{\mathit{ℬ}}^{\text{c}}+{\mathit{\omega }}_{\mathit{ℬ}}^{\text{c}}\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T\right)\right),$$(57) $${\mathit{\pi }}_{\text{m}}\equiv \frac{\partial \mathit{ℒ}}{\partial {\mathit{\omega }}_{\text{m}}}={\text{I}}_{\circ ,m}\left({\mathit{\omega }}_{\mathit{ℬ}}^{\text{m}}+f_m\left(\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T{\mathit{\omega }}_{\mathit{ℬ}}^{\text{m}}+{\mathit{\omega }}_{\mathit{ℬ}}^{\text{m}}\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T\right)\right).$$(58)

To get the position of the Cassini state, we proceed as above. The periodic oscillations δB^peri of the B-matrix are computed at the leading order and injected in the first two equations of motion whose explicit expressions read $${\dot{\mathit{\pi }}}_{\text{c}}=-g\left[{\mathit{k}}_{\text{inv}},{\mathit{\pi }}_{\text{c}}\right]-{\text{I}}_{\circ ,c}{f}_c\left[\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\left({\mathit{\omega }}_{ℬ}^{\text{c}}\right)}^2\right]-{\text{I}}_{\circ }\nu \left({\mathit{\omega }}_{ℬ}^{\text{c}}-{\mathit{\omega }}_{ℬ}^{\text{m}}\right).$$(59) $$\begin{array}{l}{\dot{\mathit{\pi }}}_{\text{m}}=-g\left[{\mathit{k}}_{\text{inv}},{\mathit{\pi }}_{\text{m}}\right]-{\text{I}}_{\circ ,c}{f}_c\left[\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,{\left({\mathit{\omega }}_{ℬ}^{\text{c}}\right)}^2\right]+{\text{I}}_{\circ }\nu \left({\mathit{\omega }}_{ℬ}^{\text{c}}-{\mathit{\omega }}_{ℬ}^{\text{m}}\right)\\ -3.{\text{I}}_{\circ }{n}^2\left[\mathit{R}{\mathit{B}}^{\prime }{\mathit{R}}^T,\mathit{S}\right],\end{array}$$(60)

with now B′ ≡ R₃(θ_p)(B + δB^peri)R₃(–θ_p). The second term in the right-hand side of ${\dot{\mathit{\pi }}}_{\text{c}}$ and ${\dot{\mathit{\pi }}}_{\text{m}}$ represents the conservative coupling between the core and the mantle generated by the pressure force of the liquid as it rotates while the third term accounts for the friction at the CMB. When adding the two equations of motion, we easily check that the total angular momentum π = π_c + π_m does satisfy the same equation of motion as in the single layer case (13), as it should be.

To solve Eqs. (59) and (60), we searched for a solution of the form R = ε_m1 i + ε_m2j + ε_m3k and ω_c = ω_p (ε_c2i – ε_c1j + k). The calculation done with TRIP shows that the in-plane obliquities ε_m1 and ε_c1 and the deviations from the Cassini plane ε_m2 and ε_c2 again are decoupled from the mantle’s rotation lag ε_m3 that is given by $${\epsilon }_{m3}=\left({\kappa }_{\zeta ,s}+{\kappa }_{\zeta ,c}+N_m{\epsilon }_{m2}\right)/\left(4{\kappa }_{22}\right).$$(61)

The other angles are solution the linearised system $$A\left(\begin{array}{c}{\epsilon }_{m1}\\ {\epsilon }_{m2}\\ {\epsilon }_{c1}\\ {\epsilon }_{c2}\end{array}\right)=-\left(\begin{array}{c}{N}_m\\ 0\\ N_c\\ 0\end{array}\right),$$(62)

where A = A_v + Α_ω + Α_ζ is a 4 × 4 matrix given by $$A_{\nu }\left[\begin{array}{cccc}{D}_{\text{m}}& F& E& -F\\ -F& D_{\text{m}}& F& E\\ E& -F& D_{\text{c}}& F\\ F& E& -F& D_{\text{c}}\end{array}\right],$$(63) $$A_{\omega }=-{\kappa }_{\omega }\left[\begin{array}{cccc}\cos \left(2\varpi -2\text{Ω}\right)& \sin \left(2\varpi -2\text{Ω}\right)& 0& 0\\ \sin \left(2\varpi -2\text{Ω}\right)& -\cos \left(2\varpi -2\text{Ω}\right)& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],$$(64) $$A_{\zeta }=\left[\begin{array}{cccc}-K_1+{\kappa }_{\zeta ,c}& K_2-{\kappa }_{\zeta ,s}& K_3& -K_4\\ -K_2+{\kappa }_{\zeta ,s}& -K_1+{\kappa }_{\zeta ,c}& K_4& K_3\\ K_3& -K_4& -K_5& K_6\\ K_4& K_3& -K_6& -K_5\end{array}\right],$$(65)

with the parameters $$E\equiv M{R}^2{C}_{20}^{\text{c}}{\omega }_p\left({\omega }_p+\dot{\text{Ω}}\cos i\right),$$(66) $$D_m\equiv \kappa +C_m{\omega }_p\dot{\text{Ω}}\cos i-M{R}^2{C}_{20}^{\text{c}}{\omega }_p\dot{\text{Ω}}\cos i-E,$$(67) $$D_c\equiv C_c{\omega }_p\dot{\text{Ω}}\cos i+M{R}^2{C}_{20}^{\text{c}}{\omega }_p\dot{\text{Ω}}\cos i-E,$$(68) $$F\equiv {\text{I}}_{\circ }\nu {\omega }_p,$$(69) $$N_m\equiv N-N_c,$$(70) $$N_c\equiv \dot{\text{Ω}}\sin i\left(C_c{\omega }_p+M{R}^2{C}_{20}^{\text{c}}{\omega }_p\right),$$(71)

where C_c and C_m are the polar moment of inertia of the core and of the mantle, respectively, and $$K_1\equiv \frac{1}{3}{k}_2{q}_r{n}^{2}M{R}^2{t}_c\left(3X_0^{-3,0}\left(e\right)+p^2{t}_c\right),$$(72) $$K_2\equiv p\zeta K_1,$$(73) $$K_3\equiv \frac{1}{6}{k}_2{q}_r{n}^{2}M{R}^2{t}_c\left(3X_0^{-3,0}\left(e\right)+p^2{t}_c\right),$$(74) $$K_4\equiv p\zeta K_3,$$(75) $$K_5\equiv \frac{1}{3}{k}_2{q}_r{\omega }_p^{2}M{R}^2{t}_c^2,$$(76) $$K_6\equiv p\zeta K_5,$$(77)

where l_c = I_∘c f_c/I_∘. The coefficients κ, κ₂₂, and κ_ω have exactly the same expression as in the single layer model (see Eqs. (17)–(19)). In this system of equations, we retained only the leading terms, disregarding those involving ${\dot{\text{Ω}}}^2$ and smaller. A keen observer might notice that the system governed by Equations (59) and (60) contains six independent equations, whereas Equations (61) and (62) together provide only five. Indeed, there is an additional equation that the system must satisfy, corresponding to the projection of ${\dot{\mathit{\pi }}}_{\text{c}}$ onto the unit vector k. This equation has no solution under the assumption that the core’s spin rate ω_c equals the mantle’s spin rate, ω_p. However, if we allow the core spin rate to differ slightly from the mantle spin rate by a small amount δω, the final equation yields $$\delta \omega \equiv {\omega }_{\text{c}}-{\omega }_p=\frac{\left(C_c+M{R}^2{C}_{20}^{\text{c}}\right){\epsilon }_{c2}{\omega }_p\dot{\text{Ω}}\sin i}{{\text{I}}_{\circ }\nu }.$$(78)

In Equation (62), we decomposed the matrix A into three matrices, A_v, Α_ω, and Α_ζ to emphasise the distinct physical processes involved. The matrix A_v encapsulates the core-mantle interaction, namely the pressure torque (parameter E) and the viscous coupling (parameter F). The matrix Α_ω describes the apsidal precession of the orbit, as observed in Peale et al. (2014) and detailed in Baland et al. (2017). The matrix Α_ζ accounts for the tidal deformation.

By neglecting apsidal precession and tidal deformation, and considering the frictionless regime (F = 0) at the CMB, the system (62) decouples into two independent subsystems involving (ε_m1, ε_c1) and (ε_m2, ε_c2), respectively. The 2 × 2 matrix with entries D_m, E, E, and D_c is generically invertible, implying that ε_m2 and ε_c2 must be zero. This signifies that, in this case, the Cassini state resides within the Cassini plane (Boué 2020). However, the inclusion of friction, apsidal precession, or tidal deformation couples (ε_m1, ε_c1) and (s_m2, ε_c2), leading to a deviation of the Cassini state from the Cassini plane.

In the following, solutions are obtained by inverting the system described in Equation (62). This analytical approximation holds valid as long as the quantities ε_m1, ε_m2, ε_m3, ε_c1, ε_c2, and δω/ω_p are much smaller than one.

5 Results

Here, we apply the model to Mercury where p = 3/2, γ₀ = 0 and σ_p = +1. We briefly consider the solid case without molten core to compare our results with those of Baland et al. (2017) and MacPherson & Dumberry (2022). Then, we rapidly switch to the two-layer model for comparison with Peale et al. (2014) and to explore a larger domain of the parameter space.

5.1 Mercury as a solid body with no fluid core

Using the parameters from Baland et al. (2017), we find the obliquity of Mercury in the Cassini plane to be ε_Ω = 2.026 arcmin (Eq. (21)) for the rigid case and ${\epsilon }_{\text{Ω}}^{k_2}=2.031$ arcmin (Eq. (46)) for the deformable case, in perfect agreement with Baland et al. (2017).

The radius of the circle described by the top of the spin axis, as the argument of the perihelion spans 180 degrees (see, e.g. Baland et al. 2017, Fig. 8), is given by our formula (23) as ε_ω = 0.898 arcsec (hereafter, we refer to this motion as the apsidal circle). The value we obtain for ε_ω is slightly larger than the 0.868 and 0.863 arcsec found by Baland et al. (2017) for the deformable and rigid cases, respectively.

However, we obtain a value closer to that of Baland et al. (2017), namely ε_ω = 0.865 arcsec, when substituting C₂₀ and C₂₂ with $C_{20}+\delta C_{20}^{\text{stat}}$ and $C_{22}+\delta C_{22}^{\text{stat}}$, respectively. As explained earlier, we consider the C₂₀ and C₂₂ values measured with MESSENGER (Mazarico et al. 2014) to include the static component of their tidal deformation. Therefore, we favour the formula without the static corrections $\delta C_{20}^{\text{stat}}$ and $\delta C_{22}^{\text{stat}}$.

Lastly, for the deviation of the apsidal circle from the Cassini plane, we obtain ε_ζ = 0.0097 arcsec, a value much smaller than its radius, in agreement with the green circle represented in Peale et al. (2014, Fig. 3), whereas Baland et al. (2017) found 0.995 arcsec, a value approximately 100 times larger. As discussed previously, we believe the extreme value reported by Baland et al. (2017) is due to incorrect projections of their equations of motion (59), which are not self-consistent because of the lack of the rotation lag (ε₃ in our notation). This omission resulted in the unit vector ŝ and its time derivative not being orthogonal.

Surprisingly, MacPherson & Dumberry (2022) reproduce the deviation angle ε_ζ from the Cassini plane, as found by Baland et al. (2017), using their Equation (16), despite working in the mantle’s equatorial plane, which by definition, is orthogonal to the figure axis. Therefore, the argument we presented against Baland et al. (2017) does not apply in this case. The origin of the discrepancy between our result and that of MacPherson & Dumberry (2022) remains unclear.

Fig. 1 Projected equilibrium spin axes on the orbital plane for varying core-mantle viscous friction, assuming a spherical core. a) Mantle spin axis coordinates. b) Core spin axis coordinates. The vertical line labelled Cassini plane in panel a) is the intersection of the Cassini plane with the orbital plane. Red dots indicate coordinates from Peale et al. (2014, Table 2). The observations are taken from Margot et al. (2012); Mazarico et al. (2014); Stark et al. (2015a); Verma & Margot (2016); Genova et al. (2019); Konopliv et al. (2020); Bertone et al. (2021).

5.2 Influence of the molten core

In Peale et al. (2014), the authors determined the equilibrium of Mercury’s spin axis through numerical integration of equations of motion similar to our system (56), but also incorporating magnetic and topographic core-mantle couplings. One of our goals is to reproduce their results analytically, eliminating the need for numerical simulations. For easier comparison, we adopt their planetary parameters but retain the orbital parameters from Baland et al. (2017).

To validate our model, we reproduced Figure 3 from Peale et al. (2014) in our Figure 1, with red dots indicating the values from their Table 2. This figure shows the projection of the equilibrium spin axis onto the orbital plane for varying core-mantle viscous friction. Our model’s agreement with their numerical solution is excellent. We also plot recent measurements of Mercury’s spin pole orientation with their associated 1σ error bars. The model can explain the observed deviation from the Cassini plane, ε_m2, provided the tilt is towards positive abscissa (X_m > 0). Among the seven measurements in Figure 1, only that of Mazarico et al. (2014) has a negative abscissa, with its error bar crossing the Cassini plane at the bottom of the figure. As noted in the introduction, the whole planet’s moment of inertia, C, can be adjusted to fit the in-plane obliquity ε_m1 without affecting the longitude libration fit, provided the mantle’s moment of inertia, C_m, remains fixed (Peale 1976). Thus, fitting the spin axis orientation primarily involves determining the core’s moment of inertia, C_c, and the viscosity, parameterised by the relaxation time τ_ν.

Nevertheless, it is important to note that our model’s leading-order computation of the planet’s elastic deformation does not account for the fluid core spin obliquity. Including this effect could cause the viscous coupling to tilt the mantle spin axis towards X_m < 0, as explained in MacPherson & Dumberry (2022).

Figure 3 of Peale et al. (2014) and our Figure 1 were obtained for a spherical core. However, Peale et al. (2014) noted that the hydrostatic shape of the core would lead to a geometric flattening of ε_f = 7.161 × 10^–5, corresponding to f_c = 1.04 in our notation. They also considered a non-hydrostatic case with twice the flattening, implying f_c = 2.08. With these values, Peale et al. (2014) concluded that the spin axes of the core and mantle should be nearly aligned and close to the Cassini plane.

To better understand this conclusion, we constructed a map of the mantle’s in-plane obliquity, ε_m1, and deviation from the Cassini plane, ε_m2, as functions of the core flattening coefficient, f_c, and viscous timescale τ_ν. The five cases considered by Peale et al. (2014) in their Figure 3 are indicated by black crosses in our Fig. 2. The hydrostatic flattening of the core places Mercury slightly outside the map to the right. In Figure 2a, we observe two plateaus: the light red area corresponds to ε_m1 ≈ 2 arcmin as observed, while the blue upper left corner shows a plateau at 1 arcmin. The ratio between these two values is exactly C/C_m, implying that only the mantle is tilted in the blue region of the parameter space (Fig. 2a). Indeed, in the upper left corner, both flattening and viscous friction are low, decoupling the core from the mantle even on the nodal precession timescale. The transition occurring at f_c ≈ 0.005 corresponds to a resonance of the precession angular speed with the nearly diurnal free wobble eigenfrequency (Ragazzo et al. 2022). The parameters chosen by Peale et al. (2014) ensure that the pressure torque is negligible (f_c ≪ 1), and the range of viscous timescales spans the transition between coupled and uncoupled core-mantle states. Figure 2b shows that this parameter space region exhibits significant deviation from the Cassini plane. As observed by Peale et al. (2014), adding the pressure torque (f_c close to one) causes the planet to behave like a solid body, reducing the deviation from the Cassini plane that is then solely induced by tidal deformation.

In Figures 2a and 2b, the white domain corresponds to extreme solution values where the linear approximation used in this study fails. The parameter space studied by Boué (2020) corresponds to the upper line of the figure (τ_ν → ∞ and a full range of f_c). Although the linear model always provides a single solution, the non-linear solution obtained in the inviscid case exhibits a bifurcation (Fig. 3). It may seem surprising that Mercury’s obliquity shows non-linear effects despite its small value. This is due to the core reaching high obliquities of the order of 90 degrees, as shown in Peale et al. (2014) and Boué (2020). In this parameter space region, it is necessary at least to treat the core’s in-plane obliquity, ε_c1, non-linearly. But as the core tilts, its angular speed significantly differs from the mantle’s, requiring the addition of a new variable ω_c to model this spin rate. And the relative angular speed, δω = ω_c – ω_p, between the core and the mantle, must also be treated non-linearly.

Fig. 2 Mantle’s in-plane obliquity εm1 (a) and deviation from the Cassini plane εm2 (b) as functions of the core flattening coefficient fc and viscous relaxation time τν. The white domains correspond to values that fall outside the colour scales. In these regions, the linear approximation is no longer valid and the equations of motions have to be solved treating the core’s in-plane obliquity and angular velocity non-linearly. Black crosses indicate the parameters considered by Peale et al. (2014), shown as red dots in our Figure 1.

6 Conclusion

In this paper, we present a new formulation of the classical model describing the rotation of Mercury, grounded in an advanced Lagrangian formalism developed by Ragazzo & Ruiz (2015) and Boué (2020). This approach enables a comprehensive 3D solution for the mantle’s orientation, including the phase in longitude at the perihelion, and in the future, will contribute to constraining its internal structure.

The final model takes into account a deformable mantle and a fluid ellipsoidal core, whose flattening imposes the conservative pressure torque. We also consider apsidal precession and dissipative effects associated with tidal interactions and friction between the core and mantle.

Firstly, we compare our results with those of Baland et al. (2017) by treating Mercury initially as a rigid body and subsequently as a deformable body. We obtain similar values for the offset in the Cassini plane, with 2.026 arcmin in the rigid case and 2.031 arcmin in the deformable case. However, we find a radius of 0.898 arcsec for the circle traced by the rotation axis during one perihelion precession, which is approximately 3% larger than the value reported by Baland et al. (2017). This minor discrepancy might be attributed to a potential double-counting of the static component of the tidal deformation effects in their approach. More significantly, our calculation of the tidally induced deviation from the Cassini plane results in a tilt that is 100 times smaller than the value reported by Baland et al. (2017). We suggest that the out-of-Cassini-plane offset may have been overestimated in their study due to an incorrect projection of the angular momentum equations onto the invariant plane.

Secondly, we obtain for the complete model the derivation of analytical equilibrium equations for Mercury’s rotation axis. We reproduce the equilibrium positions obtained numerically by Peale et al. (2014), which are compatible with the latest observations of Mercury’s rotation axis. We focus particularly on the influence of core flattening and CMB friction on both in-plane and out-of-plane offsets from the Cassini plane. We confirm the result of Margot et al. (2012), which concluded that if the core and mantle are too weakly coupled, it becomes impossible to explain the observed in-plane offset. As in Peale et al. (2014), we show that for a nearly spherical core, friction plays a dominant role, even more so than in the case of a more flattened core. This leads to a regime where the out-of-Cassini-plane offset becomes excessively large, rendering it incompatible with observations.

Finally, we identify a parameter space region where the linear approximation in Mercury’s orientation breaks down. This region is characterised by weak core-mantle friction and intermediate core flattening. In this regime, a bifurcation occurs, as described in the frictionless limit by Boué (2020), when the core angular momentum tilts beyond 30 degrees. Additionally, the core’s angular speed, forced by core-mantle friction, should significantly decrease, in contradiction with our assumption that the core and mantle angular speeds are equal.

Nevertheless, comparing our results with those of Peale et al. (2014) – which are not subject to the linear approximation – shows that a significant out-of-Cassini-plane deviation can still occur within the validity range of our model. This scenario assumes a nearly spherical core and a viscous relaxation time, driven by core-mantle friction, of the order of 10⁷ days. In this situation, tidal deformation plays a minimal role compared to the core flattening and CMB friction.

Future work is still needed to constrain Mercury’s interior structure using the upcoming data from the BepiColombo space mission. This work will focus on incorporating a solid inner core, a realistic fluid rotation model for the outer core, a more complex rheological model, and a more refined orbital motion that goes beyond the uniform precession at constant inclination and eccentricity.

Fig. 3 Mantle’s in-plane obliquity without viscous dissipation at the CMB (τν → ∞). The curve labelled ‘linear’ was computed using the present model and represents a transverse view of the top line in Fig. 2a. The non-linear curve is from Boué (2020).

Acknowledgements

We would like to thank Benoît Noyelles, whose valuable comments helped to improve the quality of the paper. C.R. is partially supported by FAPESP grant 2023/07076-4.

References

All Figures

Fig. 1 Projected equilibrium spin axes on the orbital plane for varying core-mantle viscous friction, assuming a spherical core. a) Mantle spin axis coordinates. b) Core spin axis coordinates. The vertical line labelled Cassini plane in panel a) is the intersection of the Cassini plane with the orbital plane. Red dots indicate coordinates from Peale et al. (2014, Table 2). The observations are taken from Margot et al. (2012); Mazarico et al. (2014); Stark et al. (2015a); Verma & Margot (2016); Genova et al. (2019); Konopliv et al. (2020); Bertone et al. (2021).

In the text

Fig. 2 Mantle’s in-plane obliquity εm1 (a) and deviation from the Cassini plane εm2 (b) as functions of the core flattening coefficient fc and viscous relaxation time τν. The white domains correspond to values that fall outside the colour scales. In these regions, the linear approximation is no longer valid and the equations of motions have to be solved treating the core’s in-plane obliquity and angular velocity non-linearly. Black crosses indicate the parameters considered by Peale et al. (2014), shown as red dots in our Figure 1.

In the text

	Fig. 3 Mantle’s in-plane obliquity without viscous dissipation at the CMB (τν → ∞). The curve labelled ‘linear’ was computed using the present model and represents a transverse view of the top line in Fig. 2a. The non-linear curve is from Boué (2020).
In the text