Up: Accurate spin axes and

Subsections

2 Theoretical considerations

2.1 Definition of units

In all calculations below we use the following units. The unit of time is one Julian day (1 d), where 1 d is exactly 86400 s. The sidereal day is given by 86164.09054 s (JPL). The unit of length is 1 astronomical unit (1 AU = m). Moreover the gravitational constant is chosen as G=1. In this system of units, 1 Solar mass is = .

2.2 Initial conditions and input data

Although the number of figures in numerical input data below is sometimes greater than the estimated accuracy, those are the input parameters used in our program. We present them here simply for reproducibility reasons.

All planetary masses are derived from the mass ratios listed at the JPL database (JPL). However, the Earth and the Moon are integrated separately. Their individual masses are therefore computed using the mass ratios of Sun/(Earth+Moon) and Moon/Earth. The mass ratios are listed in Table 1.

 Sun/Mercury 6023600 Sun/Venus 408523.71 Sun/(Earth+Moon) 328900.56 Sun/(Mars system) 3098708 Sun/(Jupiter system) 1047.3486 Sun/(Saturn system) 3497.898 Sun/(Uranus system) 22902.98 Sun/(Neptune system) 19412.24 Sun/(Pluto system) Moon/Earth 0.012300034 Phobos/Mars Deimos/Mars

Relative masses for Phobos and Deimos as well as the positions and velocities for all the bodies are obtained from "Horizons On-Line Ephemeris System'' (DE-406) (JPL). The reference frame is ICRF/J2000.0. All our simulations start at MJD 51600.5 (2000-Feb.-26 00:00).

The moments of inertia for the Earth are based on the JPL data (Jz-Jy)/Jz =0.0032737634, Jz/mR2 =0.3307007, R=6378.136 km, (JPL). This gives that and . For Mars, the moments of inertia are: corresponding to Jz/mR2 =0.3662, (JPL) and R=3397.2 km (Bouquillon & Souchay 1999); and from Jz/Jx =1.005741and Jz/Jy =1.005044 (Bouquillon & Souchay 1999).

The obliquity of the Earth is given by degrees at MJD 51600.5 (JPL). The present obliquity of Mars is degrees (Folkner et al. 1997). Finally, the spin angular velocities (denoted ) for Earth or Mars are listed in Sects. 2.4 and 2.7.

2.3 Treatment of the N-body problem

The N-body problem is here treated in a purely classical way. In order to keep the model simple we will not include relativistic effects, solar mass-loss or tides. These effects are studied and discussed in Appendix. Two independent engines to integrate the equations of motion are applied. First, we use the SWIFT library (Levison & Duncan 2000) which contains both the Mixed Variable Symplectic method (MVS) and the Bulirsch-Stoer method. Secondly, we have implemented a simple, but general algorithm for repeated Richardson extrapolation of Euler's method. The time step in this method is cut in half an arbitrary number of times (usually five) for each time interval and then extrapolations are repeatedly carried out in order to obtain very good positions and velocities (analogous to Romberg's method). This method is applied for the tests in the Appendix.

We have written an ANSI C program where all the above numerical methods can be used. The stability of the numerical results are tested using all methods. For appropriate time steps they all yield identical results (e.g. 15 figures in the positions are the same after thousands of years of simulation). For the final simulations, we decided to apply the MVS method because of its high speed and high stability properties (e.g. stability of the total angular momentum).

The program tests that the total mechanical energy and angular momentum are conserved at every time step. If conservation is violated the particle system is moved back to the previous time step. The integration is then attempted using half the previous time step. This procedure may then be repeated until conservation is fulfilled. The simulation then continues a certain number of steps (user defined) using the smaller time step before restoring the original time step. Another check of the accuracy is to study rapid oscillations (such as moon orbits) for peculiarities, since those are most sensitive to numerical errors. Standard modifications in order to gain numerical stability are also applied in the code, e.g. summation order. Trivial checks such as that compiler optimizations do not alter the results are also done. For every computed output quantity , we check that for two distinct simulations the numerical truncation error is negligible: . We have also investigated the sensitivity in obliquity (output) to variation in some input parameters by computing the condition number C = (relative change in output)/(relative change in input). Variations in solar mass, Earth mass, Moon mass, moments of inertia and spin angular velocity typically result in C= 10-4-10-2. This indicates that the problem is well conditioned. In present computations, we use the time step  d or lower (as low as 0.0036525 d in extreme cases). Finally, the consistency of our results is also checked by comparisons with measurements and findings of other authors, whenever possible.

2.4 The instantaneous behavior of the spin axis of a rigid body - rotation of the Earth

Instantaneous studies of spin axis precession, nutation and obliquity require a method for treating the momentary dynamics of a rotating rigid body. The approach below is general for rigid bodies and in principle exact; it can be used to study oscillations at any time scale and is still convenient for numerical simulations. We use a center of mass reference system S(x,y,z), i.e. a system that always follows the center of mass G of the body and has its axes fixed in inertial space. These axes are parallel with the external heliocentric system. is also a reference system with the origin at the center of mass, but its axes are instead fixed in the precessing body. These axes coincide with the body's principal axes so the inertia tensor is always diagonal. Furthermore, the body is considered to be a rigid (not necessarily homogeneous) oblate ellipsoid with principal moments of inertia in the general case. The actual numerical values are taken from satellite measurements (see Sect. 2.2). We shall thus ignore the phenomenon of variation of latitude and assume that the spin axis is rigidly fixed along the z'-axis of S'. In the S-system, the total external torque ( ) acting on the body w.r.t. the center of mass G is related to the angular momentum ( ) according to:

 (1)

The similarity transformation between the S- and S'-system is given by:

 (2)

where remains diagonal at all times and is the Eulerian rotation matrix in the x-convention (Goldstein 1980):

 (3)

If we define , the matrix elements of then depend on time according to:

 (4)

We thus obtain the following first order nonhomogeneous differential equation:

 (5)

By replacing dt with , we end up with the numerical expression:

 (6)

which subsequently can be improved using a few Richardson extrapolations. The explicit expression within the brackets is too long and complicated to be given here; it is fully implemented in our C-program. Initial conditions for Earth are given for simulations forward in time as follows: , where radians/d. The second factor is the ratio between a Julian day and a sidereal day (JPL). The obliquity at MJD 51600.5 is degrees (JPL). will certainly at a later stage precess away from this direction but never away from the z'-axis. is determined by . In the case of , we set and . For simulations backward in time, the initial conditions that change are: , and . also needs to be calculated at the start of the simulation (and also of course at every other time step). Note that an arbitrary vector must obey the transformation between the S and -system. In  , the basis vector is always given by (0,0,1)while in S, . We therefore have the specific relation:

 (7)

This means that given in S, we also obtain and . The last angle is simply given by in the case Jx=Jy, otherwise . Finally, the angular velocities in Eq. (4) needed to determine are approximated with and , while .

2.5 External torques acting on the planet

The exact expression for external torques from particles j acting on an ellipsoid at time t may in the body fixed system be written:

 (8)

where is pointing at an arbitrary mass element dm in the ellipsoid and is the radius vector to the particle j. However, this expression becomes far too complicated to compute, even for a spheroid. An extremely good approximation is due to MacCullagh (see Goldstein 1980, p. 228). In this approximation the potential energy of an oblate ellipsoid of mass m and a particle j is:

 (9)

where is the moment of inertia about the direction of . Jx, Jy and Jz are of course the moments of inertia w.r.t. the principal axes. The torque about the center of mass G of body i exerted by all particles j may now be written as:

 (10)

The instantaneous torques from all the other planets, the Moon and the Sun are included in the summation at every time step. In the case of Mars we are also including the moons Phobos and Deimos (see Sect. 2.10). Since is computed in , it is subsequently transformed back to S ( ) so it can be used in Eq. (6).

2.6 The summer solstice for the northern hemisphere

Suppose is the radius vector from the center of the Sun pointing at the Earth-Moon system's center of mass and is its velocity. The normal vector to the ecliptic plane is then given by . Let be the spin angular velocity of Earth. If is fulfilled we have either summer or winter solstice. It may then be checked if the angle in which case we have summer solstice. However, for simulations backward in time the condition for summer solstice instead becomes since is then replaced by . An analogous calculation can also be performed to find the northern summer solstices in the case of Mars.

2.7 Rotation of Mars

One day at Mars is slightly longer than the day for Earth, namely 24.622962 h compared to 23.93419 h (JPL). The angular velocity for Mars is therefore given by radians/d. The obliquity (the angle between the Mars equator and its instantaneous orbital plane) of Mars is given by degrees (Mars pathfinder mission, Folkner et al. 1997). Currently, the summer solstice of the northern hemisphere of Mars occurs at Dec. 16, 2000, 04h UT. Given that and are known, the equations

 = (11) = (12) = (13)

completely determine the components of . Here, ( ) is the normal vector to the orbital plane of Mars and is the heliocentric position of Mars at its summer solstice. The obtained results in radians/d are:

 2.7318578 (14) -0.3438322 (15) 5.4702752. (16)

The initial conditions for the Mars fixed system ( ) is . The Euler angle is no longer directly related to the obliquity because the Mars orbital plane has an inclination of about relative to the ecliptic. This gives that the angle between (z'-axis) and the ecliptic normal (z-axis) becomes . The last Eulerian angle is arbitrarily given by . The time derivatives are set analogously to what we did in Sect. 2.4. For simulations backwards in time the changes become as follows: , and .

2.8 The length of the summer half year

Given the semi major axis a, eccentricity e and the distance r at summer solstice for a specific year, we get from the definition of eccentric anomaly (E) that at summer solstice: . The general relation between E and the true anomaly v is:

 (17)

The specific true anomaly from perihelion to the planet at summer solstice then becomes:

 (18)

The angle from perihelion to the vernal equinox is given by and for the autumnal equinox, . The change in time dt is directly related to the change in eccentric anomaly dE through (Danby 1964):

 (19)

where n is the mean motion. Integrating between the equinoxes then gives that

 (20)

where and are the eccentric anomalies at autumnal and vernal equinoxes respectively. According to Eq. (17) we have:

 = (21) = (22)

The length of the summer half year relative to a full year is finally given by . For more discussions regarding variations of the astronomical seasons, see e.g. Berger & Loutre (1994).

2.9 Spin axis precession and orbital precession

The precession angle of the spin axis ( ) is calculated by

 (23)

where and are normal vectors to the instantaneous and present orbital plane, respectively, and analogously for and .

The orbital precession is studied by calculating the Laplace-Runge-Lenz vector which is conserved for the two-particle problem. It is directed along the radius vector to the perihelion point. In the solar system, many particle interactions (as well as the smaller general relativistic effect) cause the orbit of a planet of mass m to precess. The orbital precession can easily be studied by observing the rotating Laplace-Runge-Lenz vector which is given by:

 (24)

where and (of the planet) are the linear and angular momentum, respectively. An alternative method would be to compute the combined changes in the orbital parameters: the longitude of the ascending node , the argument of perihelion and the inclination i: d .

2.10 Torques due to Phobos and Deimos

Since Phobos and Deimos are so close to Mars (for Phobos the distance is less than 3 Mars radii), their motions are incredibly rapid compared to all other bodies in our simulation. For the time periods of interest, it is impossible to integrate these moons in the same manner as our own Moon. The objective here is therefore instead to obtain an average formula for the torque on Mars due to its moons. First a full integration will be carried out including Phobos and Deimos for a time span that is long enough to contain all relevant oscillations. The torques acting on Mars are calculated instantaneously, i.e. at every time step. We then choose circular moon orbits in the orbital plane of Mars (virtual moons) that yield exactly the same average torques as those obtained from the full integration. This treatment is similar to that given by Goldstein (1980, p. 229). Account will also be taken to the triaxiality of Mars, i.e., . The expression for I in Eq. (9) may be rewritten as:

 (25)

where the angles to a moon-position are measured relative to the body fixed system . Let us now use a coordinate system xyz for which the effective mean orbit of a virtual moon is in the xy-plane and the spin axis of Mars ( ) is in the xz-plane. The xy-plane is chosen as the instantaneous Mars orbit plane with the origin at the Mars center. The transformation between and this new xyz-system obeys: , i.e. with one obtains:

 (26)

We treat a mean virtual moon orbit as a simple circular orbit which in polar coordinates becomes: . We thus obtain:

 (27)

As usual, the moon revolution () and the rotation of Mars () are extremely rapid compared to the change in allowing a time average over an integer number of periods to be performed (the symbol refers to a time average, e.g. for a function f(t) this is given by ):

 (28)

In the same way one obtains that and . The average moment of inertia can now be written:

 (29)

The part of the averaged potential energy that has an orientational dependence is (see Eq. (9)):

 = = (30)

Finally the average torque acting on Mars due to a moon becomes:

 (31)

The torque is in the direction of , where is the normal vector to the instantaneous Mars orbit and is the spin axis of Mars. In order to apply this average formula at every time step we need to know the best possible value of r. Fortunately, we already have access to the accurate instantaneous Eq. (10). We therefore first start a simulation performing an exact integration of the Mars moons ( d). We then divide the resulting instantaneous torques by the instantaneous . By forming the average (k)

 (32)

r can be computed from . In order to get an accurate value of r, one has to be very careful in selecting the period over which the average is taken. During the simulation we found that the angle between the Mars equator and the orbital planes of Phobos and Deimos have periods of 6050 and 6050/4 years, respectively (see Figs. 1 and 2). Averages must be taken over the full periods. The resulting averages of r for Phobos and Deimos are and AU, respectively. Note that this results in much greater average torques acting on Mars than do averages taken over shorter periods. A 10 year mean treatment would instead result in and AU for Phobos and Deimos, respectively. For long simulations, it is important to use the first set of values for r.

Up: Accurate spin axes and