Latitudinal librations of Mercury with a fluid core

Icarus 203 (2009) 1–12

Contents lists available at ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

Latitudinal librations of Mercury with a fluid core Julien Dufey a,*, Benoît Noyelles a,b, Nicolas Rambaux b, Anne Lemaitre a a b

University of Namur, Department of Mathematics, 8 Rempart de la Vierge, B-5000 Namur, Belgium IMCCE, CNRS UMR 8028, Paris Observatory, Université Pierre et Marie Curie – Paris 6, USTL, 77 avenue Denfert-Rochereau, 75014 Paris, France

a r t i c l e

i n f o

Article history: Received 23 December 2008 Revised 31 March 2009 Accepted 2 April 2009 Available online 3 July 2009 Keywords: Mercury Celestial mechanics Resonances Spin-orbit Rotational dynamics

a b s t r a c t In the framework of the space missions to Mercury, an accurate model of rotation is needed. Librations around the 3:2 spin-orbit resonance as well as latitudinal librations have to be predicted with the best possible accuracy. In this paper, we use a Hamiltonian analysis and numerical integrations to study the librations of Mercury, both in longitude and latitude. Due to the proximity of the period of the free libration in longitude to the orbital period of Jupiter, the 88-day and 11.86-year contributions dominate Mercury’s libration in longitude (with the Hermean parameters chosen). The amplitude of the libration in latitude is much smaller (under 1 arcsec) and should not be detected by the space missions. Nevertheless, we point out that this amplitude could be much larger (up to several tens of arcsec) if the free period related to the libration in latitude approaches the period of the Jupiter–Saturn Great Inequality (883 years). Given the large uncertainties on the planetary parameters, this new resonant forcing on Mercury’s libration in latitude should be borne in mind. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction The rotational dynamics of Mercury is a unique case in the Solar System, because it consists of a 3:2 spin-orbit resonance, i.e. Mercury accomplishes 2 revolutions around the Sun while rotating 3 times around its spin axis. Moreover, Mercury is the target of two space missions, NASA MESSENGER, which carried out its first and second flybys of Mercury on January 14 2008 and October 6 2008, and ESA/JAXA BepiColombo, scheduled to be launched to Mercury in 2014. One of the goals of these space missions is to get information on Mercury’s internal structure, in particular its gravity field and the size of a molten core. One way to get this information is to observe Mercury’s rotation in order to detect the signature of its internal structure. In 2007, Margot et al. (2007) used radar observations to detect the signature of a fluid core in the amplitude of a forced 88-day libration in longitude. The presence of a spherical fluid core decouples the motion of the mantle from the motion of the core and external excitations such as the variation of the Sun–Mercury distance produce larger amplitudes of response. Moreover, it induces a variation of the proper frequencies of the rotation (Rambaux et al., 2007), bringing the period of the free libration in longitude close to the 11.86-year orbital period of Jupiter. Recently, three of us (Dufey et al., 2008) proposed a quasi-periodic decomposition of the expected librations in longitude, show* Corresponding author. Fax: +32 81724914. E-mail addresses: [email protected], [email protected] (J. Dufey). 0019-1035/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2009.04.005

ing that not only the 88-day libration, but also other librations due to planetary perturbations should be detected by the space missions. In particular, the contributions of Jupiter and Venus are significant. These results have been confirmed by Peale et al. (2009). In this paper, we present a two-degree of freedom theory giving not only the librations in longitude, but also in latitude, i.e. the variations of Mercury’s obliquity, with the planetary perturbations taken into account. We compute these librations applying both analytical and numerical methods. We start by giving the hypotheses of the two-degree of freedom model used in the paper. In Section 3, we describe the Hamiltonian analysis and the Lie averaging process to compute the short-periodic contributions on the librations in longitude and latitude. Section 4 presents the numerical method and the frequency analysis used to give a synthetic representation of the different variables. We state both the analytical and numerical results in Section 5 and explain the small discrepancies with an analysis of the residuals of the numerical solutions. Section 6 is dedicated to a discussion of the results. We talk about the accuracy of the determined amplitudes, the mean obliquity of Mercury, and point out a potential resonance in latitude. 2. The model The problem studied here is the problem of the rotation of Mercury with two degrees of freedom. In a previous paper (Dufey et al., 2008), we studied the planar case (one degree of freedom) with all the inclinations of the planets considered as zero. In this paper, all the orbital elements are taken into account as well as planetary perturbations acting upon them. Here are the basic hypotheses

2

J. Dufey et al. / Icarus 203 (2009) 1–12

(the units chosen are Mercury’s equatorial radius as length unit, Mercury’s mass as mass unit, and the terrestrial year as time unit):  The orbital motion of Mercury is perturbed by the other planets.  Mercury is a triaxial body (with moments of inertia A < B < C) and the gravity field is truncated after the second degree terms, i.e.

C 02 ¼ J 2 ¼

A þ B  2C 2

and C 22 ¼ C 22 ¼

BA : 4

The values of the constants are C = 0.34, C 02 ¼ 6  105 , and C 22 ¼ 1  105 (Anderson et al., 1987).  Mercury has two layers as suggested by the recent observations of Margot et al. (2007): a rigid mantle and a spherical liquid core. We neglect magnetic and viscous interaction between the mantle and the core. With this assumption, the short-term rotation of Mercury actually refers to the rotation of the mantle and, in the kinetic energy, the moments of inertia A, B, and C will be replaced by Am ; Bm ; and C m , i.e. the moments of inertia of the mantle.The ratio C m =C ¼ 0:579 and the value C 22 ¼ 1  105 used in this paper correspond to ðB  AÞ=C m ¼ 2:0319  104 which is consistent with the radar observations of Margot et al. (2007) (in their paper, the adopted value of the moment ratio is ðB  AÞ=C m ¼ ð2:03  0:12Þ  104 ).  The angular momentum axis and the third axis of inertia coincide (no wobble motion).

of the other planets (varo stands for one of the six orbital elements):

varo ¼ varH o þ fvar ðlM ; lV ; lE ; lMa ; lJ ; lS ; lU ; lN Þ;

ð2Þ

varH o

where is the Keplerian contribution of the orbital element, lM ; lV ; lE ; lMa ; lJ ; lS ; lU ; and lN are, respectively, the mean longitudes of Mercury, Venus, the barycenter Earth–Moon, Mars, Jupiter, Saturn, Uranus, and Neptune. The functions fvar are Poisson series given by the VSOP planetary theory of Simon (for a reference for the VSOP coefficients, see Fienga and Simon (2004)). In the analytical method, we will consider neither the Poisson terms of these series nor the terms depending on the longitudes of Mars, Uranus or Neptune. We check later with the numerical study that these simplifications are legitimate. Hence the series have the following form:

varo ¼ varH o þ fvar ðlM ; lV ; lE ; lJ ; lS Þ:

ð3Þ

The functions fvar are small quantities and are introduced in the Hamiltonian by series expansions. In order to keep an internal set of canonical variables, we need to make adjustments to the Hamiltonian. Since we add four variables (the longitudes of the planets) to the model and we do not consider -o and Xo as constants anymore, the new Hamiltonian is now

H¼ 3. The analytical method

l2 2L2o

_ o Ho þ _ o Go þ X þ nJ KJ þ nV KV þ nS KS þ nE KE þ -

K21 2C m

þ V G ðlo ; -o ; Xo ; eo ; ao ; io ; k1 ; k3 ; Lo ; K1 ; K3 ; lV ; lE ; lJ ; lS Þ; First, let us mention that all the analytical computations are performed using the MSNam, an algebraic manipulator. Its framework was developed by Moons (Henrard, 1986). In the body of this paper, we will only describe the main points of the analytical developments. The reader will find more details in D’Hoedt and Lemaitre (2004, 2005), Dufey et al. (2008), and in Appendices A and B.

where nJ ; nV ; nS ; and nE are the mean motions of Jupiter, Venus, Saturn, and the Earth and KJ ; KV ; KS ; and KE are the conjugated moments to lJ ; lV ; lS ; and lE . The moments Go and Ho are associated _ o are their respective precession _ o and X with -o and Xo , while rates.

3.1. The Hamiltonian and the introduction of the perturbations

3.2. Resonant angles, equilibria, and free periods

Following D’Hoedt and Lemaitre (2004) and Dufey et al. (2008), after several rotations and changes of variables, the Hamiltonian has the following form:

3.2.1. The resonant angles The first resonant angle represents the 3:2 spin-orbit resonance: r1 ¼ k1  3lo =2  -o . The second one describes the 1:1 commensurability between the orbital and rotational nodes, as stated in the third Cassini’s law (see Colombo (1966), Peale (1969) or Lemaitre et al. (2006) in the specific case of Mercury): r3 ¼ k3 þ Xo . The moments associated with r1 and r3 are, respectively, K1 and K3 , but we need to change those associated with lo ; -o ; and Xo to keep a canonical set of variables: the moment Lo will be replaced by Ko ¼ Lo þ 3K1 =2, Go by G0o ¼ Go þ K1 , and Ho by H0o ¼ Ho  K3 . The Hamiltonian is then

H¼

l2

þ 2

2L |fflffl{zfflffl}o H2B ðLo Þ

! K21 GM l3 3C 02 2 2 2 2 2  þ y Þ þ 3C ðx  y Þ ; ðx 2 2C m 2 L6o |ffl{zffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} TðK1 Þ

ð1Þ

V G ðlo ;-o ;Xo ;eo ;ao ;io ;k1 ;k3 ;Lo ;KÞ

where H2B is the Hamiltonian of the two-body problem, T is the kinetic energy, and V G is the gravitational potential expanded up to the fifth order in eccentricity, with x and y the coordinates of the unit vector ðx; y; zÞ pointing to the Sun in the frame rigidly linked to Mercury. M is the mass of the Sun, G is the gravitational constant, l ¼ GðM þ 1Þ, C m is the principal moment of inertia of the mantle, lo ; -o ; Xo ; eo ; ao ; and io are the usual orbital elements, Lo is the moment associated with the mean anomaly lo , k1 ¼ g þ h, k3 ¼ h, g is the spin angle, h is the longitude of the rotational node, K1 is the norm of the angular momentum and the moment associated with k1 , and K is the ecliptic obliquity. To be coherent with previous papers (e.g. D’Hoedt and Lemaitre, 2004), the subscript 2 refers to motions linked to the wobble, i.e. the deviation between the angular momentum axis and the polar axis. As mentioned in our hypotheses, this deviation is neglected. We now introduce indirect planetary perturbations in our model. Those perturbations do not act on the rotation itself, but on the orbit of Mercury, affecting the rotation through the 3:2 spin-orbit resonance. The orbital elements are then functions of the longitude

H¼

l2 2ðKo  32 K1 Þ2

þ nJ KJ þ nV KV þ nS KS þ nE KE

2   _ o ðH0 þ K3 Þ þ K1 _ o G0o  K1 þ X þo 2C m

þ V G ðlo ; -o ; Xo ; eo ; ao ; io ; r1 ; r3 ; Ko ; K1 ; K; lV ; lE ; lJ ; lS Þ:

ð4Þ

3.2.2. Equilibria and free periods of the averaged Hamiltonian We now average the Hamiltonian over the angular orbital elements and the planetary terms. The equilibrium of the  1 ¼ 0Þ and satisfies the Hamiltonian lies at the exact resonance ðr  3 ¼ 0Þ. The Hamiltonian used to find the third Cassini’s law ðr equilibria K1 ; K3 ; K is then

        _ o G0o  K1 H K1 ; K3 ; K ¼ H2B K1 þ T K1 þ V G ð0; 0; K1 ; KÞ þ  0  _ o H þ K3 : þX o

3

J. Dufey et al. / Icarus 203 (2009) 1–12

Then, we use several canonical changes of variables that allow us to compute the equilibria of this Hamiltonian as well as the fundamental periods of r1 and r3 :

K1 ¼ 7:70349 MM R2e years1 ; K3 ¼ 0:0577482 MM R2e years1 ;

K ¼ 7 1:197 arcmin;

T r 1 ¼ 12:055 years; and T r 3 ¼ 615:69 years; with MM the mass and Re the equatorial radius of Mercury. T r 1 corresponds to the free libration period and T r 3 to the free precession period. Details about the changes of variables and computations can be found in Appendix A. For the sake of clarity, we omit the ‘‘bars” over the variables from now on. 3.3. The generator of the Lie averaging process To compute the evolution of the variables we use a Lie averaging process (Hori, 1966; Henrard, 1970) that averages the Hamiltonian over all the angular variables. We use the Lie triangle to compute the generators of this averaging transformation up to a chosen order. Thanks to these generators, we are able to compute the evolution of any variable. The Lie triangle is

H00 H01

H10

H02 H03

H11

H20

H12

H21

H30

 with the non-averaged Hamiltonian (the input) in the first row and the averaged Hamiltonian H in the diagonal. It is obtained with the following formulas:

Hn0

¼

Hn1 1

þ



Hn1 0 ; W1



with the intermediate Hnj computed as follows:

Hnj ¼ Hn1 jþ1 þ

j    X j n1 Hji ; W iþ1 ; i i¼0

where W i is the generator of the ith floor of the Lie triangle, ( ; ) des  j ignates the Poisson bracket and is the binomial coefficient. i This process allows us to have the generators of this averaging transformation up to a chosen order. The reader will find more details about the generators in Appendix B. At this step, we use Cartesian instead of polar coordinates to avoid any singularity:

x1 ¼ r1 ; y1 ¼

K1 KH 1

x3 ¼

 1; y3 ¼

qffiffiffiffiffiffi 2K3

KH 1

qffiffiffiffiffiffi 2K3

KH 1

sinðr3 Þ; cosðr3 Þ;

with KH 1 a constant value for K1 defined in Appendix A. The formula to compute the evolution of any variable (Deprit, 1969) is

1 ; y 3 Þ þ f ðx1 ; x3 ; y1 ; y3 Þ ¼ f ðx1 ; x3 ; y

oX rder i¼1

1 ðf ðx1 ; x3 ; y1 ; y3 Þ; W i Þ; i!

ð5Þ

1 ; y 3 . with the Poisson bracket evaluated in the equilibria x1 ; x3 ; y The results of this analytical method are given in Section 5. 4. Numerical method In addition to this analytical study, we performed a numerical study in the same way as in Noyelles et al. (2008). It consists in numerically integrating the equations derived from the Hamilto-

nian (1), and then giving synthetic representations of the variables by frequency analysis. The frequency analysis algorithm that we use is named NAFF as Numerical Analysis of the Fundamental Frequencies and is based on Laskar (1993). It is improved with some refinements such as the use of different timesteps in the analysis, to detect high frequencies (see Laskar, 2003) and an iterative determination of the frequencies, taking into account the frequencies previously detected (see Champenois, 1998). The frequency analysis aims at writing a complex signal f ðtÞ under the form

f ðtÞ 

n X

ak exp

pffiffiffiffiffiffiffi  1xk t ;

ð6Þ

k¼1

where xk are real frequencies and ak complex coefficients. If the signal f ðtÞ is real, its frequency spectrum is symmetric and the complex amplitudes associated with the frequencies xk and xk are complex conjugates. This algorithm is very efficient, except when two frequencies are too close to each other. In that case the algorithm is not confident in its accuracy and stops. When the difference between two frequencies is larger than twice the frequency associated with the length of the total time interval, the determination of each fundamental frequency is not perturbed by the other ones. Although the iterative method suggested by Champenois (1998) allows to reduce this distance, some troubles still remain when the frequencies are too close to each other. The integrated equations are

dk1 oH ; ¼ oK1 dt dk3 oH ; ¼ oK3 dt

dK1 oH ; ¼ ok1 dt dK3 oH : ¼ ok3 dt

ð7Þ

The variables k0 and K0 do not appear here because the orbital motion of Mercury is determined by ephemerides and not by integration of the system. The integrations are performed with the Adams– Bashforth–Moulton 10th order predictor–corrector integrator. One of the critical points is the initial conditions. As suggested by Peale (2005), the amplitudes associated with the free librations should be negligible. Unfortunately, when choosing initial conditions even very close to the exact equilibrium, we always get a free part in the solution, that is often predominant. Several methods exist to drop this free part. For instance, Bois and Rambaux (2007) proposed to fit the mean initial conditions in order to locate the spin-orbit system at its center of libration. Another possibility is to add a damping in the equations that will reduce the amplitudes of the free librations, as done for instance by Peale et al. (2007) for the behavior of Mercury’s spin. Yseboodt and Margot (2006) explored a third way in starting from the equilibrium related to a simplified problem and in slowly switching on the planetary perturbations, to induce an adiabatic deviation of the equilibrium without apparition of free librations. In this paper we propose an iterative method, using the synthetic representation of the solutions. It consists of three steps: 1. Numerical integrations of the system over 8000 years with initial conditions ‘‘cleverly” chosen. ‘‘Cleverly” means here that we choose initial conditions as close as possible to the expected equilibrium. 2. Quasi-periodic decomposition and identification of the solutions. Since some of the solutions are diverging over the integration timespan because of very long-period contributions (especially the regression period of Mercury’s orbital node, i.e. 235,000 years), we fit and remove a polynomial to get a quasi-periodic signal. After analysis, the free part is easy to identify because the analytical study gives us a very good approximation of the

J. Dufey et al. / Icarus 203 (2009) 1–12 80

2

60

1.5

40

1

20

0.5

Node (arcsec)

Node (arcsec)

4

0 -20

0 -0.5

-40

-1

-60

-1.5

-80

-2 -4

-3

-2

-1

0

1

2

3

4

5

-4

Time (kyr)

-3

-2

-1

0

1

2

3

4

5

Time (kyr)

Fig. 1. Two numerical integrations giving the rotational node k3 : the second integration has been obtained by removing the free solution extracted from the first one by frequency analysis. This free term, with a period about 616 years, is still predominant but with an amplitude about 50 times smaller.

frequencies of the free librations. We extract the value of this free part at the time origin of the numerical simulation. The difficulty of this stage comes from the very long-period contributions. They tend to distort the quasi-periodic contributions and their period might change significantly over the integration. This alters the accuracy of the determination of the free solution at the time origin. The solution that we have found is to use the full integration to determine the polynomial to remove, but a shorter timespan to analyze the signal. The free solution can then be evaluated efficiently. 3. Removal of the free part at the time origin from the initial conditions. With that last step, we get the new initial conditions that we use in a numerical integration.

3=2lo  -o and r3 ¼ k3 þ Xo , the mean anomaly lo and the longitudes of the perihelion -o and of the ascending node Xo are linear _ o t þ /-o and functions of time, i.e. lo ðtÞ ¼ no t þ /lo , -o ðtÞ ¼ _ ot þ / _ o are constants, and _ o ; and X Xo ðtÞ ¼ X where no ; Xo /lo ; /-o ; and /Xo are the initial phases at J2000. In the following tables, we display the amplitude and the phase of each contribution: A cosðnt þ /Þ, where A is the amplitude of the contribution, n its frequency and / its phase. After applying the analytical method, our variables (we choose to display r1 ; r3 and K) are series of the following form:

Fig. 1 gives an example of use of this algorithm on the variable k3 (i.e. Mercury’s node), with sinusoidal planetary perturbations. A fourth-degree polynomial has been fitted on the evolution of k3 , then removed from it. The left panel gives the result. We can see a free sinusoidal term (v) with an amplitude of about 50 arcsec and a period around 616 years. The frequency analysis algorithm gives us an estimation of the phase and the amplitude of this term at the time origin. We removed it from the initial conditions to get the right panel, on which the free term v is still predominant, but with an amplitude of only 1.3 arcsec. A third iteration was used to analyze the evolution of this variable properly. In this third iteration, a free contribution still remains and is actually detected, but its amplitude is very small and is not considered as an obstacle in the quasi-periodic decomposition of the forced rotational motion.

where, by rotational invariance, p ¼ j þ k þ m þ n (we remind the reader that lo is the mean anomaly and not the mean longitude). We illustrate how to compute the amplitude and phase of a contribution with an example. Let us take the 11.86-year contribution: a cosðlJ  -o Þ þ b sinðlJ  -o Þ. The amplitude A is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 A ¼ a þ b and the phase is the sum of two parts: / ¼ /1 þ /2 . The first one is /1 ¼ arctanðb=aÞ and the second is the sum of the initial phases: /2 ¼ /lJ  /-o .

5. Evolution of the variables In this section we give the short-period librations of the resonant arguments r1 ; r3 and of the ecliptic obliquity K. We first present the results given by the analytical study, using series of 28 terms for the orbital ephemerides of Mercury. The series depend only on the mean longitudes of Mercury, Venus, the barycenter Earth–Moon, Jupiter and Saturn. Then, the numerical study compares the results obtained with a 28-term and a 10,000-term theory to validate the model. Finally, we discuss the accuracy of the numerical method by analyzing its residuals. 5.1. Analytical results In the analytical study, the resonant arguments are defined using the proper modes of the system. So, in the definition of r1 ¼ k1

var ¼

X

  aijkmn cos ilo þ jlV þ klE þ mlJ þ nlS  p-o

ijkmn

  þ bijkmn sin ilo þ jlV þ klE þ mlJ þ nlS  p-o ;

The short-period longitudinal librations are gathered in Table 1, while the latitudinal librations (in the ecliptic obliquity K) are in Table 3. Table 2 gives the short-period librations of the resonant angle r3 , related to the precessional motion of Mercury. Its dynamics is strongly linked to K, because they belong to the same canonical degree of freedom ðr3 ; K3 Þ, with K3 ¼ K1 ð1  cos KÞ. The longitudinal librations have already been studied in Dufey et al. (2008) and Peale et al. (2007), but with a different value of the ratio CCm . The value that we use, i.e. 0.579 is the same as in the recent numerical study of Peale et al. (2009). As they mention, the 11.862-year contribution, i.e. the orbital period of Jupiter is predominant, because of the proximity of the resonance to the free libration in longitude. Let us mention that this amplitude is strongly related to the value of the moment ratio ðB  AÞ=C m (hence to the free period). In Table 1, the value chosen is ðB  AÞ=C m ¼ 2:0319  104 . This value is consistent with the interval given by Margot et al. (2007), i.e. ðB  AÞ=C m ¼ ð2:03  0:12Þ  104 . Should we choose the inferior bound of the interval we would find an amplitude of 14.32 arcsec for the 11.86-year contribution. With the superior bound, this amplitude would be 56.94 arcsec. We also find an amplitude of 35.85 arcsec for the 88-day contribution, i.e. the orbital period of Mercury. The contributions 10–17 (except 11) have periods close to 88 or 44 days. An analytical study is very convenient to distinguish these terms

5

J. Dufey et al. / Icarus 203 (2009) 1–12

Table 1 Longitudinal librations of Mercury, obtained analytically, from the resonant argument r1 . Each contribution has the form A cosðnt þ /Þ, with A the amplitude, n the frequency and / the phase computed at J2000. No.

lo

lV

lE

lJ

lS

-o

Period

Amplitude

Ratio

Phase

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

— 1 2 2 — — 1 3 1 1 — 2 2 1 1 2 2

— — — 5 — — — — — — — — — — — — —

— — — — — — 4 — — — — — — — — — —

1  — — — 2 — — 2 2 2 1 1 1 1 2 2

— — — — 2 — — — — — 5 — — — — — —

1 — — 5 2 2 4 — 2 2 3 1 1 1 1 2 2

11.862 y 87.970 d 43.985 d 5.664 y 14.729 y 5.931 y 6.575 y 29.323 d 91.692 d 84.537 d 883.28 y 44.436 d 43.541 d 89.793 d 86.217 d 44.897 d 43.110 d

43.712 35.849 3.754 3.597 1.568 1.379 0.578 0.386 0.201 0.191 0.103 0.069 0.067 0.044 0.043 0.041 0.040

1.2193 1.0000 0.1047 0.1003 0.0437 0.0385 0.0161 0.0108 0.0056 0.0053 0.0029 0.0019 0.0019 0.0012 0.0012 0.0011 0.0011

164.46° 84.80° 79.59° 166.85° 85.96° 125.91° 25.57° 105.62° 58.69° 48.28° 153.53° 25.51° 4.71° 17.94° 7.17° 63.89° 43.07°

as as as as as as as as as as as as as as as as as

Table 2 Librations of r3 , related to the precessional motion. The first line, written in italic, represents the amplitude of the librations due to the precessional motion. Since it is a longperiod perturbation, a model studying both the core and the mantle of Mercury would be more appropriate than ours. As a consequence, this line is just an indication that longperiod effects might alter the determination of some contributions. Each contribution has the form A cosðnt þ /Þ, with A the amplitude, n the frequency and / the phase computed at J2000. No.

lo

lV

lE

lJ

lS

-o

Xo

Period

Amplitude

Ratio

Phase

1 2 3 4 5 6 7 8 9 11 10 12 13

— 1 2 — 3 4 — 2 1 — 2 — 1

— — — — — — — — — — 5 — —

— — — — — — — — — — — — 4

— — — 2 — — 2 — — 1 — — —

— — — 5 — — — — — — — 2 —

2 — 2 3 2 2 2 — 2 1 5 2 4

2 — 2 — 2 2 — — 2 — — — —

63315 y 87.970 43.985 883.280 29.323 21.992 5.931 43.985 87.970 11.862 5.664 14.729 6.575

3.74 143.06 67.31 65.44 45.02 14.71 14.64 12.91 8.20 7.83 6.00 4.43 1.93

26.128 1.0000 0.4705 0.4574 0.3147 0.1028 0.1023 0.0902 0.0573 0.0547 0.0420 0.0309 0.0135

56.84° 84.80° 67.25° 15.01° 107.55° 77.66° 117.10° 100.41° 62.04° 148.63° 100.24° 157.27° 66.00°

from the first 88 and 44-day ones, but it might not be so easy for radar observations (Margot et al., 2007) or observations by the space missions MESSENGER and BepiColombo. The librations related to the latitudinal motion (Tables 2 and 3) are dominated by a very long-period contribution (63315 years) corresponding to twice the argument of the perihelion ð2xo Þ. It results in an adiabatic oscillation of Mercury’s spin axis, that is beyond the scope of this paper because we study only the shortperiod oscillations. Moreover, as mentioned by Peale (2006), the

d d y d d y d d y y y y

as mas mas mas mas mas mas mas mas mas mas mas mas

adiabatic deviation of Mercury’s spin axis is followed by the mantle and by the molten core, so our decoupling hypotheses are not valid for low frequencies. As a consequence, the amplitude that we indicate for this contribution should not be considered as reliable, but only as an indication that long-period effects might alter the determination of some contributions. However, we are confident in the amplitudes associated with the other contributions. Peale (2006) has already suggested that the short period oscillations of Mercury’s obliquity would be smaller than 1 arcsec. Table 3

Table 3 Latitudinal librations, i.e. librations around the equilibrium of the ecliptic obliquity K ¼ 7:01995 . Each contribution has the form A cosðnt þ /Þ, with A the amplitude, n the frequency and / the phase computed at J2000. No.

lo

lV

lE

lJ

lS

-o

Xo

Period

Amplitude

Ratio

Phase

1 2 3 4 5 6 7 8 9 10 11 12 13

— 2 3 — 1 4 — 2 — 1 2 — 1

— — — — — — — — — — 5 — —

— — — — — — — — — — — — 4

— v — 2 — — 2 — 1 — — — —

— — — 5 — — — — — — — 2 —

2 2 2 3 — 2 2 — 1 2 5 2 4

2 2 2 — — 2 — — — 2 — — —

63315 y 43.985 29.323 883.280 87.970 21.992 5.931 43.985 11.862 87.970 5.664 14.729 6.575

457.230 8.213 5.551 4.676 3.599 1.798 1.624 1.483 1.259 0.991 0.523 0.477 0.185

127.0400 2.2820 1.5304 1.2991 1.0000 0.4996 0.4511 0.4120 0.3499 0.2754 0.1452 0.1325 0.0515

33.16° 22.75° 162.45° 134.22° 174.79° 12.34° 144.13° 169.59° 108.13° 27.96° 24.59° 100.25° 171.72°

d d y d d y d y d y y y

mas mas mas mas mas mas mas mas mas mas mas mas mas

6

J. Dufey et al. / Icarus 203 (2009) 1–12

confirms this result by giving amplitudes smaller than 10 mas. We can notice that the 44-day and 29-day contributions, i.e. twice and thrice Mercury’s orbital frequency, are larger than the 88-day one. More interesting is the 883-year contribution, due to the Jupiter– Saturn Great Inequality. This contribution is due to the proximity of the Jupiter–Saturn system to the 5:2 orbital resonance. It has a significant dynamical effect on the Solar System bodies, on the planets themselves (Laplace, 1785), on the Asteroidal Belt (Ferraz-Mello et al., 1997; Henrard, 1997) and also on the irregular satellites of the outer planets (C´uk and Burns, 2004). In addition to this effect, one can notice contributions of Jupiter (5.931 and 11.862 years), Saturn (14.723 years), Venus (5.664 years) and the Earth (6.575 years).

that are not included in the series, a new Model 1 is elaborated to include these contributions. The frequencies are estimated independently for each model. The longitudinal librations (Table 4) are obtained after substraction of a polynomial from the spin angle k1 , to remove the mean spin rate of Mercury (39.13 rad y1) and the contributions due to the very long periods. Then we obtain a signal that we assume to be quasi-periodic, and, using the method described in Section 4, we analyze it. As explained earlier, the frequency analysis algorithm determines first the contributions associated with the highest frequency, then the next one, and so on, until the detected frequency is very close to one already determined. The detection of such a frequency is generally due to a lack of accuracy in the determination of the frequencies. Here, the slow variations of g o induce time variations of the detected frequencies, hence a loss of accuracy in the determination of frequencies. That is the reason why the determination stops before the detection of the smaller contributions. We still find a very good agreement between numerical and analytical results. When we compare Models 1, 1b and 2, we can notice a dramatical change of the amplitude associated with lJ , i.e. the 11.862-year contribution. The closeness to the resonance makes this amplitude very sensitive to the parameters of the system, as shown by Peale et al. (2009). This explains a lack of accuracy in its determination. The librations associated with the latitudinal motion (Tables 5 and 6) have been obtained in a similar way. They tend to confirm the analytical results, especially for the presence of the 883-year contribution. However, the comparison between the Models 1 and 1b shows significant differences for the amplitudes associated with the lo and 2lo contributions (respectively 87.97 and 43.98 days). This can be explained by the presence, in the analytical results (Tables 2 and 3), of significant contributions in lo ; lo þ 2g o ; 2lo and 2lo þ 2g o that the frequency analysis is unable to distinguish. So, the variation of the amplitude that is numerically derived is due to the slow evolution of g o . In Model 1b, the 883-year contribution might be quite difficult to determinate because the time length of the interval of study is not large compared to 883 years. Moreover, this contribution has to be distinguished from a 616-year one, i.e. a residual of the free precession. Finally, a very good agreement between the Models 1 and 2 can be seen. It validates the choice of the 28 terms used in the Model 1 and the analytical model.

5.2. Numerical results We give here our numerical results, obtained after removal of the free librations. Three analyses have been made, named as Models 1, 1b, and 2.  In Model 1, the frequency analyses have been performed over 262,144 points with a timestep of 10.5192 days, the starting point being JD 990559.4639 (i.e. 4000 years before J2000), and the time interval 7500 years. The ephemerides are modeled with the series of 28 terms used in the analytical study.  In Model 1b, the frequency analyses are performed over 65,536 points with a timestep of 21.0384 days, the starting point being JD 2358265.8476 (i.e. 255 years before J2000), and the time interval 3800 years. Since the time interval is shorter than in Model 1, we expect to detect less frequencies. The ephemerides are the same as in Model 1.  In Model 2, the frequency analyses are the same as in Model 1, but the ephemerides are modelized with about 10,000 terms. These terms are all quasi-periodic terms present in the VSOP ephemerides. Model 1 aims at comparing a numerical and an analytical study of the same model. Model 1b, with a different starting date, shows the time-dependency of the amplitudes. Since the slow variations of g o cannot be isolated by frequency analysis, the amplitudes found are likely to be time-dependent. The aim of Model 2 is to validate the choice of the 28 terms used in the analytical study. If the frequency analyses of the results of Model 2 show contributions

Table 4 Longitudinal librations obtained numerically from

r1 , defined as k1 from which a polynomial has been removed.

No.

lo

lV

lE

lJ

lS

-o

Period

1 2 3 4 5 6

— 1 2 2 — —

— — — 5 — —

— — — — — —

1 — — — — 2

— — — — 2 —

1 — — 5 2 2

11.862 87.970 43.985 5.664 14.729 5.931

Table 5 Numerical computation of the librations of

y d d y y y

Model 1

Model 2

Model 1b

42.443 35.834 3.772 3.565 — —

41.786 35.864 3.772 3.593 1.572 1.377

47.557 as 35.671 as – — — —

as as as as

as as as as as as

r3 , associated with the precessional motion.

No.

lo

lV

lE

lJ

lS

-o

Period

1 2 3 4 5 6 7

1 2 — 3 4 — 2

— — — — — — 5

— — — — — — —

— — 2 — — 2 —

— — 5 — — — —

— — 3 — — 2 5

87.970 43.985 883.280 29.323 21.992 5.931 5.664

d d y d y y y

Model 1

Model 2

Model 1b

143.59 68.65 65.30 47.57 15.41 14.95 5.71

143.63 68.65 65.30 47.57 15.41 14.65 5.83

146.67 71.85 65.93 48.82 — — —

mas mas mas mas mas mas mas

mas mas mas mas mas mas mas

mas mas mas mas

7

J. Dufey et al. / Icarus 203 (2009) 1–12 Table 6 Latitudinal librations, obtained after numerical studies of the ecliptic obliquity K. No.

lo

lV

lE

lJ

lS

-o

Period

1 2 3 4 5 6

2 1 3 — — —

— — — — — —

— — — — — —

— — — 2 2 1

— — — 5 — —

— — — 3 2 1

43.985 87.970 29.323 883.280 5.931 11.862

5.3. Analysis of the residuals The analytical results give more frequencies than the numerical ones. To find the missing contributions, we compute the Fourier spectra of the residuals of our numerical solutions. We obtain the residuals by removing the series given by Tables 4–6 from the solutions obtained after the numerical integrations. First, the Fourier spectra indicate significant amplitudes for the contributions with periods close to harmonics of the mean motion, i.e. 87.97 days, 43.98 days, 29.32 days, and so on. Unfortunately, it is impossible to isolate all these contributions to identify them. So, the Fourier spectra at this range of periods does not give much information (see Fig. 2). On the contrary, plotting them on a longer timescale allows to estimate the contribution of several planetary perturbations that the frequency analysis failed to detect. The Fourier spectra of

d d d y y y

Model 1

Model 2

Model 1b

9.617 6.428 6.104 4.706 1.659 1.261

9.617 6.428 6.104 4.706 1.631 1.262

9.466 5.727 6.078 4.997 — —

mas mas mas mas mas mas

mas mas mas mas mas mas

mas mas mas mas

Models 1/1b are given in Fig. 3 and Table 7. Here, we detect for instance the 6.6-year (due to the Earth) and 3.9-year ð3lJ Þ contributions that are detected analytically but not numerically. The Fourier spectra from Model 2 are given in Fig. 4 and Table 8. They present quite the same visual aspect as Model 1 for the range of periods 1-18 years, so we plot them only on 18–90 years, showing the periodic contributions that are neglected in the Model 1. A lot of information appears in these spectra. For instance, we can see two contributions that are close to 2lJ (5.931 years), i.e. 5lS (5.891 years) and 4lJ  5lS (5.971 years). The proximity of these two contributions to 2lJ could explain the difference obtained between Models 1 and 2 for the nodal motion (r3, Table 5). We can also notice the anecdotical contribution of Mars in 1.881 year. More interesting are the 84 and 42-year contributions. In the table they are labelled as due to Uranus, but we should not forget

Fig. 2. Two examples of residuals of the frequency analysis, for the longitudinal librations (left), and for the latitudinal librations (right) of Model 1.

0.6

0.9 0.8

0.5 0.7 0.4

milli-arcsec

arcsec

0.6 0.5 0.4 0.3

0.3

0.2

0.2 0.1 0.1 0

0 2

4

6

8

10

Period (yr)

12

14

16

18

2

4

6

8

10

12

14

16

18

Period (yr)

Fig. 3. Fourier spectra of the residuals of the frequency analysis, for the longitudinal librations (left), and for the latitudinal librations (right) of Model 1. Here, only the contributions with periods between 1 and 18 years are plotted. The Fourier spectra for Model 2 have the same visual aspect.

8

J. Dufey et al. / Icarus 203 (2009) 1–12

Table 7 Fourier spectra of the residuals of the solutions given by the 28 terms in the Model 1. lo

lV

lE

lJ

lS

-o

Period

1 1 — 2 1 — —

2 3 — 5 — — —

— — — — 4 — —

— — 3 — — 1 —

— — — — — — 2

2 3 3 5 4 1 2

1.110 1.380 3.954 5.664 6.575 11.862 14.729

0.09

y y y y y y y

r1

r3

K

63 mas 5 mas 56 mas — 578 mas 0.7 as —

0.57 mas 0.88 mas 0.98 mas 1:6  102 mas 1.9 mas 5.9 mas 4.1 mas

7:4  102 mas 7:5  102 mas 0.11 mas 0.51 mas 0.18 mas 2:5  102 mas 0.47 mas

0.25

0.08 0.2

0.07

milli-arcsec

arcsec

0.06 0.05 0.04

0.15

0.1

0.03 0.02

0.05

0.01 0

0 20

30

40

50

60

70

20

80

30

40

50

60

70

80

Period (yr)

Period (yr)

Fig. 4. Fourier spectra of the residuals of the frequency analysis, for the longitudinal librations (left), and for the latitudinal librations, both in the Model 2 (i.e. 10,000-term ephemerides). Here, only the contributions with periods between 18 and 90 years is plotted. The Fourier spectra for Models 1/1b do not give significant information in this range.

Table 8 Fourier spectra of the residuals of the solutions given by the complete VSOP orbital ephemerides (10,000 terms). lo

lV

lE

lm

lJ

lS

lu

-o

Period

1 1 — — — 2 — — — 1 — — — — — — — — — —

2 3 — — — 5 — — — — — — — — — — — — — –

— — — — — — — — — 4 — — — — — – — — — -

— — 1 — — — — — — — — — — — — — — — — —

— — — 4 3 — — 2 4 — — — 1 — 1 — 2 — 1 —

— — — — — — 5 — 5 — 4 3 — 2 1 1 4 — 2 —

— — — — — — — — — — — — — — — — — 2 — 1

2 3 1 4 3 5 5 2 1 4 4 3 1 2 — 1 2 2 1 1

1.110 1.380 1.881 2.965 3.954 5.664 5.891 5.931 5.971 6.575 7.364 9.819 11.862 14.729 19.859 29.457 30.473 42.010 60.947 84.020

that the system Uranus–Neptune is close to a 2:1 mean-motion resonance. These two contributions cannot be distinguished in our Fourier spectra. 6. Discussion 6.1. Accuracy of the determined amplitudes A good determination of the amplitudes of the periodic contributions is crucial to get information of the internal structure of

y y y y y y y y y y y y y y y y y y y y

r1

r3

K

59 mas 12 mas — — 60 mas — 9.2 mas 8.2 mas 9.1 mas 568 mas — 136 mas 380 mas — 21 mas 85 mas — 52 mas — 16 mas

0.64 mas 0.95 mas 6:7  103 mas 5:9  102 mas 1.01 mas 1:46  102 mas 8:4  102 mas 0.15 mas 8:9  102 mas 1.8 mas 2:2  102 mas 0.31 mas 5.9 mas 4.2 mas 5:3  102 mas 0.30 mas 5:8  102 mas 0.63 mas 0.11 mas 0.75 mas

7:2  102 mas 8:6  102 mas 6:8  104 mas 6:4  103 mas 0.11 mas 0.49 mas 1:15  102 mas 9:0  103 mas 8:6  103 mas 0.18 mas 2:5  103 mas 3:6  102 mas 2:75  102 mas 0.46 mas 6:6  103 mas 0.20 mas — 6:9  102 mas 4:6  102 mas 8:3  102 mas

Mercury. For instance, Margot et al. (2007) use the 88-day longitudinal libration to constrain the size of Mercury’s core. Our study highlights some sources of error that might occur in the planetological interpretation of observed librations. We stress in particular two librations for which the interpretation might be altered:  87.97 days (the orbital period of Mercury). Because of its orbital eccentricity, longitudinal librations of the same period appear (Comstock and Bills, 2003). Unfortunately, there are several contributions with periods close to 88 days. Five of them

9

J. Dufey et al. / Icarus 203 (2009) 1–12 0.8

0.4

Obliquity (as)

Obliquity (as)

0.6

0.2 0 -0.2 -0.4 -0.6 -0.8 0

500

1000

1500

2000

25 20 15 10 5 0 -5 -10 -15 -20 -25 0

500

Time (years)

1000

1500

2000

Time (years)

Fig. 5. Left: variation of the obliquity around its equilibrium for T r3 ¼ 615, 853, and 880 years. Right: same evolution, with T r3 ¼ 880 and 883 years.

ðlo  2lJ and lo  lJ Þ are listed in Table 1, with periods of 84.5, 91.7, 86.2 and 89.8 days. The accumulated amplitude of these contributions is around 0.5 arcsec. Hence, in addition to the observational error, an other one of maximum 0.5 arcsec (when all the contributions are in phase) should be taken into account in inverting the observed amplitude.  5.93 years (half of the orbital period of Jupiter). Like the first two ones, this contribution is of high interest in the orbital and rotational dynamics of Mercury. The Fourier spectra related to the node and obliquity of Mercury show two contributions with periods very close to 5.93 years, i.e. 5lS (5.89 years) and 4lJ  5lS (5.97 years) (see Table 8). These perturbations are, respectively, due to the orbital motion of Saturn, and to Jupiter–Saturn interactions. Their estimated amplitudes, relatively to the one associated with 2lJ (see Tables 2 and 3) are about 1%. These additional contributions should be taken into account when using the observed amplitude associated with 5.93 years. In addition to these problems of contributions too close in the frequency space, the rotational dynamics of Mercury depends on several very long-period contributions, especially in r3 and obliquity. These contributions, due to the precessional motion of the node and the perihelion of Mercury, are the signature of relativistic effects on Mercury’s orbit (Einstein, 1915). As stated in Section 5, these perturbations induce a time-dependency on the amplitudes that can alter not only a numerical prediction (see Tables 4–6) but also an observation.

gravity coefficients of Mercury are not well known. The J 2 is known with an uncertainty of 33% and C 22 presents an uncertainty of 50%. In this range of values, the free period T r3 could very well be 883 years or close to it. In this case, the rotational motion of Mercury could be strongly modified. Let us compute the evolutions of r3 and the obliquity with different values of J 2 and C 22 , resulting in a change of the free period T r3 . Fig. 5 shows the variation of the obliquity around its equilibrium, with different values of T r3 obtained with the analytical method. The amplitude of the libration in 883.2 years can be very large (more than 20 arcsec for T r3 ¼ 883 years). The evolution of the resonant angle r3 presents a similar behavior and is not represented here. Fig. 6 shows the contribution of the Great Inequality on the ecliptic obliquity K in function of the free period T r3 . Here again, we see that when T r3 comes close to 883.2 years, the amplitude of the libration due to the Great Inequality tends to be very large. Even though the ratio ðB  AÞ=C m is better known with the radar measurements of Margot et al. (2007), there is still some uncertainty on the gravity coefficients J 2 and C 22 . Another resonant forcing is then possible when T r3 is close to 883 years (within the range of uncertainty of J 2 and C 22 ), and produces variations up to several tens of arcsec in the latitude. Let us note that due to small ellipticity at the core–mantle boundary, the long-period librations of the core tend to follow the mantle and in this case, the spherical core assumption is not valid anymore (gyrostatic rigidity, Poincaré, 1910). 7. Conclusion

6.2. About the mean obliquity of Mercury

6.3. A potential resonance in latitude In our model, the core is spherical and the resonant period T r3 is 615.69 years. As highlighted in Section 5, the oscillation in 883.2 years, resulting from the Great Inequality between Jupiter and Saturn, is the fourth largest contribution in the obliquity. The proximity between this forced period and the free period T r3 increases significantly the amplitude of the libration (inversely proportional to the difference between these two frequencies). Moreover the

Through analytical and numerical approaches, we have studied the longitudinal and latitudinal librations of Mercury. We take the

Amplitude of the great inequality libration (as)

Recently, Margot et al. (2007) confirmed by radar analysis that Mercury occupies the Cassini state 1. In particular, the normal to its orbit, the spin orientation and the normal to its Laplace Plane are coplanar. A mean obliquity with respect to the orbital plane of 2.11 ± 0.1 arcmin has been measured. The value for the equilibrium of the ecliptic obliquity K found in Section 3.2.2 does not match the equilibrium measured. The reason is that we consider a model with a spherical fluid core, valid on the short-term and K is the equilibrium of this model. Our study cannot be used to characterize the Cassini state and should be seen as complementary to the long-term ones of Peale (2006), Yseboodt and Margot (2006), Bois and Rambaux (2007), Rambaux et al. (2007) and D’Hoedt and Lemaitre (2008).

12 10 8 6 4 2 0 800

820

840

860

880

900

920

940

Period of σ3 (years) Fig. 6. Importance of the Great Inequality contribution on the ecliptic obliquity K in function of the period of the free libration in latitude.

10

J. Dufey et al. / Icarus 203 (2009) 1–12

updated values of the physical parameters of Mercury (Margot et al., 2007) and for a moment ratio of ðB  AÞ=C m ¼ 2:0319  104 , our results for the longitudinal librations are in good agreement with those of Peale et al. (2009). Since one of the free periods of the rotation (12.06 years) is now close to Jupiter’s orbital period, the libration in longitude is dominated by Jupiter’s contribution (11.86-year period), with an amplitude around 45 arcsec and by the 88-day period contribution, with an amplitude around 36 arcsec. Our theory also allows us to predict the latitudinal librations of Mercury. We expect these librations to be well under 1 arcsec and should therefore not be detected by the space missions. However, we explore the possibility of a resonant forcing on this libration. The main planetary perturbation on the libration in latitude is that of Jupiter–Saturn Great Inequality, with a period of 883 years. Similarly to the libration in longitude, if the other free period of the rotation comes close to the Great Inequality period, the amplitude of the libration in latitude could grow up to several tens of arcseconds. A good way to study this resonance would be to model the response of Mercury to mid-period perturbations, i.e. around 1000 years. It would also be interesting to model the inertial coupling between the core and the mantle for such long periods. Finally, to have a complete model of rotation, the wobble motion (i.e. the angle between the angular momentum and the third axis of inertia) should also be included in the future.

The authors were supported by the contract Prodex C90253 ‘‘ROMEO” from BELSPO.

sffiffiffiffiffiffiffiffiffi 2K3 x3 ¼ sinðr3 Þ; KH1 sffiffiffiffiffiffiffiffiffi 2 K3 y3 ¼ cosðr3 Þ: KH1

ðA2Þ

These new variables are canonical with multiplier 1=KH 1. One could wonder why we choose those particular y1 ; x3 ; and y3 pffiffiffiffiffiffiffiffiffi 0 0 instead of the more immediate K1 ¼ KH 2K3 cosðr3 Þ, 1 þ y1 , x3 ¼ pffiffiffiffiffiffiffiffiffi 0 and y3 ¼ 2K3 cosðr3 Þ. The reason is that K3 does not appear as so in the Hamiltonian, but through cos K and sin K. The relations between K3 and K are

K3 cos K ¼ 1  ; K1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin K ¼ 1  cos2 K : We have

¼1

 K3 K3 K3  ¼ 1  H 1  y1 þ y21  y31 þ    ¼1 H K1 K1 ð1 þ y1 Þ K1  x23 þ y23  1  y1 þ y21  y31 þ    : 2

A similar development using y01 ; x03 ; and y03 would yield the following:

Appendix A. Changes of variables In this appendix, we show how to compute the equilibria and the free periods of the Hamiltonian averaged over the angular orbital elements and the planetary terms. A.1. The averaged Hamiltonian H and the value KH 1 Averaging the Hamiltonian (4) over lo ; -o ; Xo and the longitudes of the other planets gives the following averaged Hamiltonian:1

K2 _ o ðG0o  K1 Þ þ 1 þ V G ðr1 ; r3 ; K1 ; KÞ þ 2 3 2C m 2ðKo  2 K1 Þ

l2

0 1 !2 !3 0 0 x23 þ y23 @ y01 y01 y01 cos K ¼ 1  1 Hþ2 6 þ    A: 2KH K1 KH1 KH1 1

where the constants have been evaluated at the values given in the VSOP theory, and the constant terms dropped. Since H does not depend on the mean anomaly, the associated moment Ko is a constant. To compute the equilibria of this averaged Hamiltonian we must use Cartesian coordinates. First we need to define a special value KH 1. Assuming that Mercury is a perfect sphere with a rigid mantle and a spherical liquid core, the Hamiltonian is simply the sum of the two-body contribution and the kinetic rotational energy: K21 l2 H ¼  2L 2 þ 2C . The 3:2 spin-orbit resonance and Hamilton equam o tions give us a special value for the angular momentum K1 :

ðA1Þ

where no is the mean motion of Mercury.

1 All the variables in this averaged Hamiltonian should be barred to note that they are mean variables, but for the sake of clarity, we omit them.

ðA3Þ

The expression is much simpler with the variables that we choose. Furthermore, these variables are dimensionless, which is always more careful. As in Dufey et al. (2008), expanding the two-body problem Hamiltonian H2B and the kinetic energy T around KH 1 gives

H2B þ T ¼ 

_ o ðH0 þ K3 Þ; þX o

3 KH1 ¼ no C m ; 2

To compute the equilibria of the averaged Hamiltonian, we use the following change of variables: ðr1 ; K1 ; r3 ; K3 Þ ! ðx1 ; y1 ; x3 ; y3 Þ, with r1 ¼ x1 (series expansions of cosðr1 Þ and sinðr1 Þ), K1 ¼ KH 1 ð1 þ y1 Þ and

cos K ¼ 1 

Acknowledgments

H¼

A.2. The change of variables to Cartesian coordinates

2

l2 2ðKo  32 KH 1Þ

þ 2

KH1 y21 ; 2C m

ðA4Þ

with the linear term in y1 vanishing. We also omit the first term since it is a constant in the averaged Hamiltonian H. After these expansions and changes of variables, the averaged Hamiltonian can be written

1 H¼ H K1

order X

! ai;j;k;l xi1 yj1 xk3 yl3

;

ðA5Þ

iþjþkþl¼0

where ‘‘order” is the order of the series expansion. A.3. The equilibria and the free periods We already know that at the equilibrium (exact resonances),

rH1 ¼ rH3 ¼ 0, hence xH1 ¼ xH3 ¼ 0 and Eq. (A5) becomes 1 Heq ðy1 ; y3 Þ ¼ H K1

order X

! aj;l yj1 yl3

:

ðA6Þ

jþl¼0

With a simple iterative process, it is now possible to find the equilibria for y1 and y3 :

11

J. Dufey et al. / Icarus 203 (2009) 1–12

with the intermediate Hnj computed as follows:

9 yH 1 ¼ 4:1937474  10 ;

yH 3 ¼ 0:12244468; giving 2

2

xH þ yH 3  K ¼ arccos 1  3H  2K1 1 þ yH 1 H

! ¼ 7:01995

¼ 7 1:197 arcmin:

ðA7Þ

We perform a translation to put these equilibria at the center of our coordinates. To find the fundamental periods of r1 and r3 we take the quadratic averaged Hamiltonian:

Hð2Þ ðx1 ; y1 ; x3 ; y3 Þ ¼ 0:34466  102 x21 þ 0:40274  103 x1 x3 þ 0:50814  102 x23 þ 19:566y21  0:61460  103 y1 y3 þ 0:49220  102 y23 : We need to dissociate the two degrees of freedom, i.e. perform another canonical transformation that removes the mixed terms in x1 x3 and y1 y3 . For this, we use the untangling transformation (Henrard and Lemaitre, 2005) and after a last change to action-angle variables

(

pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2U sinðuÞ; y1 ¼ 2U cosðuÞ; pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi x3 ¼ 2V sinðvÞ; y3 ¼ 2V cosðvÞ; x1 ¼

ðA8Þ

the quadratic Hamiltonian is now

Hð2Þ ð; ; U; VÞ ¼ nU U þ nV V;

ðA9Þ

with nU and nV the free frequencies of the degrees of freedom 1 and 3. The corresponding free periods are:

T r1 ¼ 2p=nU ¼ 12:055 years; T r3 ¼ 2p=nV ¼ 615:69 years: Appendix B. Computing the generators of the transformation of the Lie averaging process The Lie triangle is

H00 H01 H02 H03

H10 H11

H20

H12

H21

H30

with the non-averaged Hamiltonian (the input) in the first column, P i H ¼ ni¼0 H0i i! (the small parameter  being implicit or explicit) and P i the averaged one in the diagonal: H ¼ ni¼0 Hi0 i! . 0 The choice of the different Hi is decisive. For example, if we tried to compute the generators of the averaging transformation without making the changes of variables of the previous section, we would have extremely small denominators, hence, no convergence of the process. Our choice for H00 is

l2

H00 ¼   2 þ nu U þ nv V þ nV KV þ nE KE þ nJ KJ 2 Ko  32 KH 1 þ nS KS : H01 ,

ðB1Þ

we put the periodic terms of the Hamiltonian. P i To obtain H ¼ ni¼0 Hi0 i! , we use the following homological equations:

  Hn0 ¼ Hn1 þ Hn1 1 0 ; W1

where W i is the generator of the ith floor of the Lie triangle. Let us show how it works for the first and second order. The first order homological equation is

H10  H01 ¼

oH00 oW 1 oH00 oW 1 oH0 oW 1  þ0 0 þ0 olo oKo oKo olo oU ou |{z} 0



oH00 oW 1 oH0 oW 1 oH0 oW 1 þ0 0 þ0 þ0 0 oV ov oKV olV oKE olE

oH00 oW 1 oH0 oW 1 oW 1 oW 1 þ0 0 ¼ no  nu oKJ olJ oKS olS olo ou oW 1 oW 1 oW 1 oW 1 oW 1  nE  nJ  nS :  nv  nV ov olV olE olJ olS 

We choose H01 to have only periodic terms. Consequently, H10 ¼ 0 since it is the average of H01 , and H10  H01 ¼ H01 . From the last equation, it is now possible to compute the generator of the first order W 1 by simple integrations over the angular variables. Thanks to the generator W 1 we are now able to compute the second floor of the Lie triangle. First we have to compute the inter    1 1 ðH01 ; W 1 Þ þ ðH00 ; W 2 Þ. However, we mediate H11 ¼ H02 þ 0 1 have two unknowns in this formula: H11 and W 2 . To solve this, we introduce this equation in the second homological equation H20 ¼ H11 þ ðH10 ; W 1 Þ and we have

      1   1 H20 ¼ H02 þ H01 ; W 1 þ H10 ; W 1 þ H00 ; W 2 : 0 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} e H 02

ðB2Þ

e 0 , the previous equation allows us Since it is possible to compute H 2 to compute the generator of this second order W 2 . The other orders are computed in the same way. References



In

j    X j Hn1 ji ; W 1þi ; i i¼0

Hnj ¼ Hn1 jþ1 þ

Anderson, J.D., Colombo, G., Esposito, P.B., Lau, E.L., Trager, G.B., 1987. The mass, gravity field, and ephemeris of Mercury. Icarus 71, 337–349. Bois, E., Rambaux, N., 2007. On the oscillations in Mercury’s obliquity. Icarus 192, 308–317. Champenois, S., 1998. Dynamique de la résonance entre Mimas et Téthys, premier et troisième satellites de Saturne, Ph.D. Thesis, Observatoire de Paris. Colombo, G., 1966. Cassini’s second and third laws. Astron. J. 71, 891–896. Comstock, R.L., Bills, B.G., 2003. A Solar System survey of forced librations in longitude. J. Geophys. Res. 108 (E09), 5100. C´uk, M., Burns, J., 2004. On the secular behavior of irregular satellites. Astron. J. 128, 2518–2541. Deprit, A., 1969. Canonical transformations depending on a small parameter. Celest. Mech. 1, 12–30. D’Hoedt, S., Lemaitre, A., 2004. The spin-orbit resonant rotation of Mercury: A two degree of freedom Hamiltonian model. Cel. Mech. Dyn. Astr. 89, 267–283. D’Hoedt, S., Lemaitre, A., 2005. The spin-orbit resonance of Mercury: A Hamiltonian approach. In: Kurtz, D.W., Bromage, G.E. (Eds.), Transits of Venus: New Views of the Solar System and Galaxy, Proceedings of IAU Colloquium 196, held 7–11 June, 2004 in Preston, UK, pp. 263–270. D’Hoedt, S., Lemaitre, A., 2008. Planetary long periodic terms in Mercury’s rotation: A two-dimensional adiabatic approach. Cel. Mech. Dyn. Astr. 101, 127–139. Dufey, J., Lemaitre, A., Rambaux, N., 2008. Planetary perturbations on Mercury’s libration in longitude. Cel. Mech. Dyn. Astr. 101, 141–157. Einstein, A., 1997. The foundation of the general theory of relativity. In: Klein, M.J., Kox, A.J., Schulmann, R., (Eds.), The Collected papers of Albert Einstein. The Berlin years: Writings 1914–1917, Doc. 30, Vol. 6. Princeton University Press, pp. 146–200. Ferraz-Mello, S., Nesvorny´, D., Michtchenko, T.A., 1997. On the lack of asteroids in the Hecuba Gap. Cel. Mech. Dyn. Astr. 69, 171–185. Fienga, A., Simon, J.-L., 2004. Analytical and numerical studies of asteroid perturbations on solar system planet dynamics. Astron. Astrophys. 429, 361– 367.

12

J. Dufey et al. / Icarus 203 (2009) 1–12

Henrard, J., 1970. On a perturbation theory using Lie transforms. Celest. Mech. 3, 107–120. Henrard, J., 1986. Algebraic manipulation on computers for lunar and planetary theories. In: Kovalevsky, J., Brumberg, V. (Eds.), Proceedings of the IAU Symposium, vol. 114, Reidel, pp. 59–62. Henrard, J., 1997. The effect of the Great Inequality on the Hecuba Gap. Cel. Mech. Dyn. Astr. 69, 187–198. Henrard, J., Lemaitre, A., 2005. The untangling transformation. Astron. J. 130, 2415– 2417. Hori, G., 1966. Theory of general perturbations with unspecified canonical variables, Publications of the Astronomical Society of Japan, vol. 18, No. 4, pp. 287–296. Laplace P.S., 1785. Théorie de Jupiter et de Saturne, Mémoires de l’Académie royale des Sciences de Paris. Laskar, J., 1993. Frequency analysis of a dynamical system. Cel. Mech. Dyn. Astr. 56, 191–196. Laskar, J., 2003. Frequency map analysis and quasiperiodic decompositions. In: Proceedings of Porquerolles School, math/0305364. Lemaitre, A., D’Hoedt, S., Rambaux, N., 2006. The 3:2 spin-orbit resonant motion of Mercury. Cel. Mech. Dyn. Astr. 95, 213–224.

Margot, J.-L., Peale, S.J., Jurgens, R.F., Slade, M.A., Holin, I.V., 2007. Large longitude libration of Mercury reveals a molten core. Science 316, 710. Noyelles, B., Lemaıtre, A., Vienne, A., 2008. Titan’s rotation: A three-dimensional theory. Astron. Astrophys. 478, 959–970. Peale, S.J., 1969. Generalized Cassini’s laws. Astron. J. 74, 483–489. Peale, S.J., 2005. The free precession and libration of Mercury. Icarus 178, 4–18. Peale, S.J., 2006. The proximity of Mercury’s spin to Cassini state 1 from adiabatic invariance. Icarus 181, 338–347. Peale, S.J., Yseboodt, M., Margot, J.-L., 2007. Long-period forcing of Mercury’s libration in longitude. Icarus 187, 365–373. Peale, S.J., Margot, J.-L., Yseboodt, M., 2009. Resonant forcing of Mercury’s libration in longitude. Icarus 199, 1–8. Poincaré, H., 1910. Sur la précession des corps déformables. Bull. Astron. 27, 321– 356. Rambaux, N., Lemaitre, A., D’Hoedt, S., 2007. Coupled rotational motion of Mercury. Astron. Astrophys. 470, 741–747. Rambaux, N., van Hoolst, T., Dehant, V., Bois, E., 2007. Inertial core–mantle coupling and libration of Mercury. Astron. Astrophys. 468, 711–719. Yseboodt, M., Margot, J.-L., 2006. Evolution of Mercury’s obliquity. Icarus 181, 327– 337.