The Role of Secondary Resonances in the Evolution of the Mimas–Tethys System

Icarus 140, 106–121 (1999) Article ID icar.1999.6115, available online at http://www.idealibrary.com on

The Role of Secondary Resonances in the Evolution of the Mimas–Tethys System Sylvain Champenois and Alain Vienne Laboratoire d’Astronomie, Bureau des Longitudes, 1 Impasse de l’Observatoire, F-59000 Lille, France E-mail: [email protected] Received October 20, 1997; revised January 15, 1999

We have numerically investigated the role of the 200-year-long period, recently discovered, on the past tidal evolution of the Mimas–Tethys system through the present inclination-type resonance. We show that it has a deciding effect on the descriptions of the resonance motion, as considered up to now. Several modes of evolution are found possible according to the value allowed for the present eccentricity of Tethys (which is very badly known). As a result, Mimas’s inclination before capture may have been higher (up to 0.7◦ ) or lower (down to 0.03◦ ) than the value previously considered (0.42◦ ). Also, Tethys’s eccentricity on capture may have been quite higher (≈0.008 versus 0). Moreover, calculation of the probability of capture is found to be more complicated. We develop a method which takes chaos into account, and high probabilities of capture are found for quite small values of Tethys’s eccentricity at capture time (e.g., a probability of 0.76 for an eccentricity of 8 × 10−4 ). °c 1999 Academic Press Key Words: chaos; secondary resonances; tides; capture probability.

slowly oscillates around zero with a great amplitude (95◦ with a period of 70 years). This is the greatest libration amplitude amongst the major satellites in the Solar System. Since the argument 2φ is present in terms factored by i 1 i 3 in the expansion of the perturbative potential, we shall call this resonance the i 1 i 3 resonance. Since it is the center of our study, we shall sometimes also call it the “major resonance.” The caracteristic of this resonance is the oscillation of the direction of the conjunction of Mimas and Tethys with Saturn about the mid-point of the ascending nodes of the two satellites. In fact, there are six primary resonances at the 2 : 4 mean motion commensurability, whose arguments are as follows: 2λ1 − 4λ3 + 2Ä1 2φ = 2λ1 − 4λ3 + Ä1 + Ä3 2λ1 − 4λ3 + 2Ä3 λ1 − 2λ3 + $3 2λ1 − 4λ3 + $1 + $3

1. INTRODUCTION

λ1 − 2λ3 + $1 .

Mimas and Tethys are two satellites orbiting Saturn with respective periods of about 1 and 2 days. Their orbital elements are referred to the mass center and equatorial plane of Saturn and to the ascending node of this plane on the J2000 ecliptic. We use the following notations for the osculating elements: a, n, e, i, γ , $, Ä, and λ are the orbital semimajor axis, mean motion, eccentricity, inclination, sine of semiinclination, longitude of periapse, longitude of ascending node, and mean longitude, respectively. As we shall in fact consider the six main inner satellites of Saturn, the subscripts 1, 2, 3, 4, 5 and 6 will refer, throughout the paper, to Mimas, Enceladus, Tethys, Dione, Rhea, and Titan, in that order. We set αi j = ai /a j and (m i )i=1,6 as the satellites’ masses in units of Saturn’s mass. The Mimas–Tethys pair is strongly connected by a 2 : 4 inclination-type resonance that is unique in the Solar System: the argument 2φ = 2λ1 − 4λ3 + Ä1 + Ä3

The tides cause the satellites to evolve on converging orbits, so that the system encounters the resonances in the order listed above and may possibly be captured in one of them. As a matter of fact, the internal satellite (Mimas) moves away from the planet faster than the external satellite (Tethys), thus increasing the ratio α = a1 /a3 of their semimajor axes, and decreasing the value of 2n 1 − 4n 3 , which equals in turn the average values ˙ 1 −Ä ˙ 3 , −2Ä ˙ 3 , −2$ ˙ 1 , −Ä ˙ 3 , −$ ˙1−$ ˙ 3 , and −2$ ˙ 1 . It can −2Ä be deduced that before getting captured in the present i 1 i 3 resonance, Mimas and Tethys first encountered the i 12 resonance, from which they escaped. These various resonances do not overlap, because of the large oblateness of Saturn and the small masses of its satellites, contrary to what occurs for the corresponding resonances in the system of the satellites of Uranus (Dermott et al. 1988). Besides the resonance, Mimas and Tethys are mainly disturbed by this large oblateness of Saturn. That is why their orbits

106 0019-1035/99 $30.00 c 1999 by Academic Press Copyright ° All rights of reproduction in any form reserved.

107

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

were often described by slowly precessing ellipses whose semimajor axes, eccentricities and inclinations were assumed to be constant (the precession periods are on the order of a year). Moreover, Tethys was supposed to move on a circular orbit. In this case, the resonance brings only periodic terms in the mean ˙ 6φ, ˙ and 10φ). ˙ longitudes (with frequencies near to 2φ, At present, the ephemerides of this pair of satellites, as well as those of the six other major satellites of Saturn, are computed with the theory TASS1.6 (Vienne and Duriez 1995, Duriez and Vienne 1997). Contrary to the previous theories, TASS is constructed in a dynamical consistent way, in which the satellites are considered all together; its parameters are all independent; they are explicitly the initial conditions, the masses of the satellites, and the mass and the oblateness coefficients of Saturn. The internal precision of TASS is a few tens of kilometers over one century. The root-mean-square residuals of TASS over one century of Earth-based observations reach 00 0,12 for the best data sets, until 00 0,015 for the mutual events. The solutions in TASS for the motion of Mimas and Tethys appear very different from all the previous theories. There are more numerous terms in the series which are far from negligible. Some of them have periods greater than 100 years and up to 800 years (Vienne and Duriez 1992). So, compared with the fundamental frequencies of the system (from 1 and 2 days for the mean motion to 1 and 5 years for the precession motions), these periods are very long. Since the discovery of these long-period terms, we have conducted some partial studies about them. In Vienne et al. (1992), we showed that the acceleration in the longitude of Mimas computed by Kozai (1957) and Dourneau (1987) cannot be accounted for by tidal effects as they have suggested. In fact, we showed that this “acceleration” is a numerical artifact due to the fact that these authors did not take account in their series of the very long periods we have discovered. In Vienne et al. (1996), we began the study of the secular resonance and their eventual chaotic features: a variation of about 10−5 rad year−1 over 200,000 years of the frequency libration has been detected for initial conditions located near the separatrix of the phase ˙ which is an indicator of chaotic motion (Laskar plane (2φ, 2φ), 1990). Among the long period terms we have discovered, there is one which is particularly interesting, because of its period (200.15 years). We denote its frequency by σ˙ . In Champenois and Vienne (1999, henceforth CV) we investigate the role it has played in the tidal evolution of the system using a perturbed-pendulum model, showing that it may have brought chaos to the system, which changes the values of e3 and i 1 before capture as derived by Allan (1969). We also noticed the theoretical possibility of a capture into the 1/1 secondary resonance between σ˙ and the libration frequency of the major resonance, which makes possible a series of captures and escapes of the major resonance over an extended period, instead of a single capture. We confirm these theoretical considerations in the present paper thanks to numerical integrations, and we extend them to a

capture in the 1/2 secondary resonance. We also find the variation of the frequency σ˙ to be less than 8% since capture. Moreover, we show that the computation of the probability of capture, as proposed by Sinclair (1972), is no more valid because of the presence of the chaotic layer. A method is proposed in this paper, and we find high probabilities for quite small values of the eccentricity of Tethys. The frequency σ˙ is found to issue from the following arguments: ψ + 2φ = 3φ + σ ψ = φ+σ ψ − 2φ = −φ + σ, where σ =

Ä1 3Ä3 − + $3 . 2 2

Their degree in eccentricity and inclination is 3, which is only one more than that of the term with argument 2φ corresponding to the major resonance (the arguments defined with the next harmonics of 2φ would be at least of degree 5 and can therefore be neglected). Let us notice that the terms of these arguments are factored by i 1 , i 3 , and e3 in the expansion of the perturbative potential. Hence the eccentricity e1 of Mimas does not intervene, which explains that in spite of its nonnegligible value (about 0.02), this quantity has little effect on the long-term evolution of the system. In effect, the periods of the arguments which contain the pericenter of Mimas are too small (about 1 year against 70 for the libration period). The quantities i 1 and i 3 for the present time are fairly well known from TASS1.6 (i 1 = 1.62 ± 0.02◦ , i 3 = 1.093 ± 0.003◦ ). On the other hand, the present value e30 of e3 is badly known. The TASS1.6 theory sets it at e30 = 0.000235 but with such uncertainties that it can amount to 0.001. This is a very small value, which certainly explains why the previous studies (Allan 1969, Sinclair 1972, and Dermott et al. 1988) would assume Tethys moving on a circular orbit. We show that so small a value nevertheless exerts a deciding influence on the dynamics of the Mimas–Tethys system. Further˙ 1−Ä ˙ 3 +$ ˙3 more, it remains that the combination n 1 − 2n 3 + Ä is very near to zero whatever the value of e30 may be. Our study finally differs from the similar ones conducted about the tidal evolution of the uranian satellites (Tittemore and Wisdom 1988, 1989, Malhotra and Dermot 1990) in that we consider third-order resonances. In the particular case of the saturnian system, we may average over the second-order resonances (except the Mimas–Tethys one), since they are isolated (but this is no more true to third order (CV)). On the contrary, the above authors cannot apply the averaging principle to second-order resonances in the uranian system because of the overlapping. Besides, we differ from the authors who studied the evolution of the jovian satellite system (Malhotra 1991, Showman and Malhotra 1997) by the type of resonance we are considering:

108

CHAMPENOIS AND VIENNE

we deal with a two-body resonance (perturbed by the other saturnian satellites), whereas these authors deal with the three-body Laplace resonance (perturbed by Callisto). First of all, we recall in Section 2 the major caracteristics of Allan’s (1969) and Sinclair’s (1972) analyses about the dynamics of the Mimas–Tethys system under the hypothesis of circular orbits. The former, integrating backward in time, finds the following values for the satellites inclinations before capture in the present resonance: i 1 = 0.42◦ , i 3 = 1.05◦ . The latter finds the probability of capture to be 0.04. In Section 3 we investigate the role that tides may have played on the evolution of the system. In addition to the tidal effects on the mean motions, considered by Allan, we also take into account the effects of tides in Tethys, which tend to circularize the orbit of this satellite. We also present in this section the physical model we use for our integrations. In Section 4, we integrate backward in time our equations for different values of e30 and notice three major possible scenarii: in the first one, the value of i 1 before the system was captured in the present i 1 i 3 resonance ranges from 0.45◦ to 0.70◦ according to e30 , and quite large values of e3 are found possible at capture time (up to 0.008). In the second one, we found i 1 ≈ 0.18◦ and e3 = 0.008 before capture. In that case, the system got trapped in the 1/2 secondary resonance between the libration period of the i 1 i 3 resonance and the 200-year-long period. A third scenario was found to be possible, in which there was a capture in the 1/1 secondary resonance, which means that i 1 ≈ 0.03◦ and e3 ≈ 0 before capture. Furthermore, we found the frequency σ˙ to be only slowly varying between capture time and present time (about 8%). All of this shows that Allan’s hypothesis (1969) of a circular orbit of Tethys needs to be reconsidered, because it does not allow for chaos and secondary resonances. In Section 5, we present surfaces of section of the Mimas– Tethys system and explain thanks to them the results obtained in Section 4. In Section 6, we generalize the formula used by Sinclair (1972) to evaluate the probability of capture in the major resonance. It is shown from numerical integrations that the value 0.04 is no more correct as soon as e3 6= 0. The probability is found to be very sensitive to the value of e30 . However, it appears to be much higher than the value found by Sinclair.

2. DYNAMICS OF THE RESONANCE FOR A CIRCULAR ORBIT OF TETHYS

Allan (1969) and Sinclair (1972) investigated the dynamics of Mimas and Tethys on capture into the major i 1 i 3 resonance. Our study being an extension to theirs, we recall in this section the major features of their work. 2.1. Allan (1969) Allan investigated the tidal evolution of the system with the two following hypotheses:

—The evolution is mainly dominated by the terms containing the argument 2φ in the expansion of the right-hand sides of Lagrange’s planetary equations (for Allan estimates its period to be much longer than that of all the other arguments). —Mimas and Tethys are assumed to move in circular orbits. Using a Legendre polynomial expansion for the disturbing function (that is not adequate for mutual attraction between satellites), he gets the secular variations of the metric elements through µ ¶ dn 1 dn 1 ∂U ∗ = 6m 3 α13 n 21 + dt ∂φ dt T µ ¶ ∗ dn 3 dn 3 ∂U = −12m 1 n 23 + dt ∂φ dt T (1) ¤ ∗ di 1 m 3 α13 n 1 £ 2 ∂U 1 + 4γ1 =− dt sin i 1 ∂φ ∗ £ ¤ di 3 ∂U m 1n3 1 − 8γ32 =− , dt sin i 3 ∂φ where U ∗ stands for the resonant part of the disturbing function (which contains terms whose arguments are multiples of φ). The suffix T indicates the variation due to tidal effects. Neglecting 2 d 2 Ä1 + ddtÄ2 3 , he gets the following equation for φ: dt 2 £ ¤ φ¨ = 6n 21 α13 m 3 + 24m 1 n 23 f 0 (α13 )γ1 γ3 sin 2φ + F,

(2)

where µ F=

dn 1 dt

µ

¶ T

dn 3 −2 dt

¶ .

(3)

T

In resonance, the mean value (taken on a libration cycle) of the ¨ is zero. As a consequence, Eq. (2) second derivative of φ, hφi, yields f 0 (α13 )γ1 γ3 hsin 2φi =

−F . 6n 21 α13 m 3 + 24n 21 m 1

(4)

Hence he gets the variation of metric elements on capture into the i 1 i 3 resonance: " m2 À À ¿ ¿ µ ¶ 1 dn 1 1 dn 3 27n 1 ae 5 13/3 251 + = =− 2 m1 n 1 dt n 3 dt 4Qs a1 + 22/3 À ¿ 1 + 4γ12 m 3 n 1 α13 F 1 di 1 =− 2 sin i 1 dt 6n 1 α13 m 3 + 24n 23 m 1 sin2 i 1 À ¿ 1 di 3 m1n3 F 1 − 8γ32 =− 2 . sin i 3 dt 6n 1 α13 m 3 + 24n 23 m 1 sin2 i 3

m 23 45 m3 42/3

#

(5)

F being negative (were it not for the resonance, tidal effects would cause Mimas to move nearer to Tethys), we can see that the inclinations increase because of tides in the planet.

109

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

This increase of the inclinations results in a decrease in the amplitude which allows Allan to estimate the time of capture, by computing the time when the amplitude was 180◦ . He thus showed that the system entered libration as a was 99.217% of its present value and that the values of i 1 and i 3 before capture were 0.42◦ and 1.05◦ , respectively. In the future, the inclinations will go on increasing, while the amplitude goes on decreasing, but these variations will be lower and lower. 2.2. Sinclair (1972) Under the same hypotheses as Allan, Sinclair numerically estimates the probability of capture in the major resonance by integrating the following equations: dn 1 = 6m 3 α13 n 21 f 0 (α13 )γ1 γ3 sin(2φ) + dt dn 3 = −12m 1 n 23 f 0 (α13 )γ1 γ3 sin(2φ) + dt

µ µ

dn 1 dt dn 3 dt

¶ ¶

T

T

di 1 1 (6) = n 1 α13 m 3 γ3 f 0 (α13 ) sin(2φ) dt 2 di 3 1 = n 3 m 1 γ1 f 0 (α13 ) sin(2φ) dt 2 µ ¶ 1 γ3 γ1 dφ = n 1 − 2n 3 + n 1 α13 m 3 + n 3 m 1 f 0 (α13 ) cos(2φ). dt 8 γ1 γ3 The inclinations he considered before capture are the ones evaluated by Allan (i 1 = 0.42◦ and i 3 = 1.05◦ ). The method he used is recalled in Section 6.1. He finds a probability of 0.04 and shows that it is independent of the value taken for Saturn’s tidal dissipation function Qs . 3. MODEL FOR ECCENTRIC ORBITS

3.1. Tidal Effects The attraction of the tidal bulge raised on Saturn by one of its satellites results in a loss of angular momentum by the planet for the benefit of the satellite. This causes the angular rotation velocity of the planet to slow down and, subsequently, the orbit of the satellite to expand (if it is above synchronous height). The secular variation induced by this physical process on the mean motion n of a satellite with semimajor axis a and relative mass m is given by (Kaula 1964) µ ¶ µ ¶ 9k2 n 2 m ae 5 dn =− + O(e2 ), (7) dt T 2Qs a where k2 and ae are Saturn’s second-order Love number and equatorial radius, respectively. We shall neglect the O(e2 ) term as we are mainly concerned with the time interval since capture in the major resonance. As a matter of fact, e1 and e3 remain weak throughout this period (see below). Gavrilov and Zharkov (1977) give k2 = 0.341 and for Qs the lower bound: Qs = 14,000. We

shall assume here that Qs is independent of both the amplitude and the frequency of the tides raised on Saturn by its satellites, so we assume it to be constant throughout Mimas’s and Tethys’s orbital evolution through the present i 1 i 3 resonance. Its exact value is unknown, but it is in direct ratio to the time (see, e.g., Dermott et al. 1988). We shall therefore express our results in Qs years rather than in years, as Malhotra did, for example, in her study of the Laplace resonance (Malhotra 1991). We can easily deduce from Eq. (7) that tidal effects on the semimajor axis are in inverse ratio to the power 11/2 of the semimajor axis and in direct ratio to the mass of the satellite. Tethys is about 17 times more massive than Mimas, but the closeness of the latter to the planet wins, so the orbit of Mimas is expanding faster than that of Tethys when no resonance ties the satellites together. Permanent capture into resonance has thus been possible as the ratio of the mean motions came near to 2. The resonance now preserves this ratio by transferring from the inner satellite (Mimas) to the outer one (Tethys) part of the angular momentum coming from the reduced speed of Saturn’s spin. The expansion of Mimas’s orbit has thus been reduced, whereas the expansion of Tethys’s has been sped up. Tidal effects in the planet also act to increase the eccentricities of the satellites if their motion is prograde (which is the case of Mimas and Tethys). However, these effects, in Saturn’s system of satellites, is quite negligible with respect to tides raised on the satellites by the planet (Peale et al. 1980). The latter tends, on the contrary, to decrease the eccentricity, so that the satellite eventually moves on a circular orbit. The decrease in the eccentricity e of a satellite is given by µ ¶ e de (8) =− dt T τ with, for a satellite in synchronous spin state (which is the case of Mimas and Tethys), τ=

38a 2 r Q, 63n 3 ρ R 4

(9)

where a, n, R, ρ, r , and Q are the semimajor axis, mean motion, radius, density, rigidity, and tidal dissipation function of the satellite, respectively (Dermott et al. 1988). τ is the eccentricity damping time scale. In the calculation of τ , all the parameters are well known, except for r and Q. Taking Mimas and Tethys for icy satellites, Dermott et al. give for Mimas τ1 = 3 × 106 Q1 years and for Tethys τ3 = 106 Q3 years. The Q of a satellite depends on the material, the tidal frequency, and the temperature. According to Nakamura and Abe (1977), the Q of an icy satellite for frequencies of a few hertz at temperatures ranging from 100 to 160 K, which we can expect in the small, icy satellites of Saturn, is about 100 to 200. However, if the satellites were colder, then Q3 may be greater. Taking for Tethys Q3 = 100 , we get τ3 = 108 years. For Mimas, however, a value of Q1 = 100 is not in accordance with its present high eccentricity (e1 ≈ 0.02 (Table I)). So, unless

110

CHAMPENOIS AND VIENNE

TABLE I Present Conditions for the Mean Motions, Inclinations, and Eccentricities of Mimas and Tethys (TASS1.6, Vienne and Duriez 1995)

Mimas Tethys

n (rad/yr)

i (deg)

e

2422.44 1213.17

1.62 1.093

0.0194 <0.001

it has been enhanced by previous eccentricity resonance(s), now disrupted, one may reasonably deduce that Q1 is in fact much bigger than 100 and that the minimum value of the eccentricity damping time scale is of the order of 109 years. Peale et al. (1980) gives the minimum value τ1 ≈ 3 × 109 with the assumption that, 4.6 × 109 years ago, e1 was of the order of 0.1. As the age of the Mimas–Tethys resonance is less than 2.4 × 108 years for Qs = 14,000 (CV), we shall neglect the tidal damping of e1 . Anyway, its role in the long-term evolution of the Mimas–Tethys system is limited to the J2 terms of the equations (see Appendix). The above value of τ3 means that the eccentricity of Tethys was, on capture into the major resonance, about 10 times higher than its present value (CV). However, if capture into the major resonance took place sooner (Qs > 14,000), then the ratio could be higher than 10. Of course, a higher value of Q3 would result in a weaker value of the eccentricity on capture, but the rigidity of Tethys may also be considerably less than that of ice due to the existence of a liquid core (Peale et al. 1980) and to the fact that cometary impacts may have reduced its cohesive strength (Dermott et al. 1988). This would then decrease τ3 very much and lead to a bigger value for the past eccentricity of Tethys.

sify the angle-dependent arguments present in the expansion of the interaction potential into three classes: the first ones have periods of the order of the day (e.g., the period of revolution of Mimas around Saturn), the second ones have periods of the order of the year (e.g., the precession period of Mimas’s node), and the third group have periods of the order of the century: these are, to third degree in eccentricity–inclination, the resonant argument 2φ (with a period of 70 years at present and more in the past (CV)) and the arguments 3φ + σ, φ + σ , and −φ + σ (with a period of about 200 years at present; see Introduction). What we call short-period terms are the first two classes. We neglect them on account of the averaging principle and keep only the four arguments 2φ, 3φ + σ, φ + σ , and −φ + σ among all the angle-dependent arguments. The secular terms of the type γ1 γ3 cos(Ä1 − Ä3 ) are thus averaged out, but we keep the angleindependent secular terms (e.g., e12 or γ12 ) up to third degree in eccentricity–inclination. Lowest-order and lowest-degree secular perturbations from Enceladus, Dione, Rhea, and Titan were also included, because of their nonnegligible contribution to the precession motions of Mimas and Tethys (again, only angle-independent secular terms have been retained, the other being averaged out). At last, we took account of both tidal dissipation in Saturn and in Tethys, with Qs = 14,000 and Q3 = 100 (see Section 3.1). However, some integrations were done using the value Q3 = 385 (we are precise when presenting the results (Section 4.2)). The system we use for our numerical integrations is described in detail in the Appendix. The pertinent physical data are given in Tables I and II. 4. NUMERICAL INTEGRATIONS

3.2. The Model The physical model was taken to be Mimas and Tethys on eccentric orbits inclined on the equatorial plane of Saturn. Saturn’s gravitational momenta are essential as they provide the main contribution to the orbital precession rates, and we had to take account in the equations of the lowest-degree oblateness terms J2 , J4 , J22 , J6 , J2 J4 , and J23 (for J2 , the next degree in eccentricity–inclination also had to be taken into account). The mutual interaction potential of Mimas and Tethys was truncated to third-degree in eccentricity–inclination and averaged over the short-period terms. There are in fact three orders of magnitude in the periods present in the Mimas–Tethys system, which clas-

4.1. The Method Our equations (see Appendix) were integrated backwards in time. The initial conditions are taken from TASS1.6 for J2000. i 1 , i 3 , and φ are fixed to the following values: i 10 = 1.62◦ , i 30 = 1.093◦ , φ 0 = 29.965◦ . Different values were chosen for n 1 and n 3 around the possible values given by TASS, so that the present libration amplitude is equal to the observed value 95 ± 1◦ . The value chosen for σ 0 is not essential, since σ makes a complete turn in about 200 years, so changing its value only means considering, for example, J1900 or J2100 instead of J2000. Therefore different values were also taken for this angle.

TABLE II Parameters of the Saturnian System up to Titan 106 m 1

106 m 2

106 m 3

106 m 4

106 m 5

106 m 6

M¯ /MS

ae

J2

J4

J6

0.0634

0.15

1.060

1.963

4.32

236.638

3498.790

60330

0.016298

−0.000915

0.000095

Note. From TASS1.6 for the satellites masses m 1 , m 3 , m 4 , and m 5 (in units of Saturn’s mass) and the oblateness coefficient J6 . The coefficients J2 and J4 , the mass m 6 of Titan, the mass MS of Saturn (given by M¯ /MS ) and its equatorial radius ae (in kilometers) are from Campbell and Anderson (1989). The mass m 2 of Enceladus is an averaged value issued from the determinations of TASS1.5 (0.107 × 10−6 ), TASS1.6 (0.069 × 10−6 ), Harper and Taylor (1993; 0.213 × 10−6 ), and Dourneau (1987; 0.206 × 10−6 ).

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

In fact, the main change in the initial conditions between the various runs concerns e30 , which was set in turn to e30 = 0, e30 = 0.000235 (TASS1.6 value), and e30 = 0.001 (the maximum possible value). The exact values for n 1 , n 3 , σ , and e3 are given explicitly in Section 4.2 before we present the results obtained from each run. We chose to accelerate somewhat the tidal evolution of the system by dividing the values of Qs and Q3 that we have chosen by the “acceleration factor” g = 5 (for some runs, we chose g = 3.85; we give in Section 4.2 the value chosen for each run). It is of course important that the ratio (Qs /Q3 ) be kept the same. These choices correspond to at most 16-CPU-h runs on our workstation (166 Mhz). Of course, accelerating tidal effects may affect the dynamical evolution of the system: if the tidal dissipation rates are near or beyond an adiabatic limit (which depends on the resonance considered), then the probabilities of capture into a given resonance may not be the same (see Malhotra and Dermott 1990). However, our calculations show that, as far as the primary i 1 i 3 resonance is concerned, the tidal dissipation rates, enhanced with the value g = 5, are still 1000 times below the adiabatic limit. Concerning secondary resonances, however, the rates with g = 5 are close to the adiabatic limit. However, we are not interested in this paper in capture probabilities into secondary resonances, but only in the different possibilities for the past evolution of the Mimas–Tethys system and in capture probabilities into the primary resonance affected by chaos (in Section 6). Furthermore, capture probabilities into secondary resonances decrease as tidal effects are increased, when the adiabatic limit is approached (Malhotra and Dermott 1990). Hence if we notice a capture into a secondary resonance for enhanced tidal dissipation rates, then we know that the possibility exists, a fortiori, for lower tidal dissipation rates. The values chosen for g allowed us to carry out our runs over at most 60 × 106 years, with an Adams 10 routine using a predictor–corrector scheme and a step length of 0.1 year. The integrator was tested over this period by comparing the results obtained with two different runs (each one with e3 ≡ 0 and no tidal effects to avoid a possible chaotic behavior), using for the second one-half the step length of the first one. After 60 × 106 years of integration, the difference on the fast variable φ is less than 1φ = 5◦ at the end of this period (that is, after more than 600,000 libration periods of 2φ). The predictor–corrector methods we have used therefore appear to be quite reliable over the period we have investigated.

111

i 1 i 3 resonance. However, we notice that capture in this resonance takes place all the later as e30 is high, as we can see on Fig. 1a, where the variation of the libration amplitude after capture can be seen for three different runs: those obtained

4.2. Results We found three possible scenarii for the past tidal evolution of the Mimas–Tethys system: The first one is characterized by n 01 = 2422.5072 rad year−1 , n 03 = 1213.1732 rad year−1 , σ 0 = 180.64◦ , Q3 = 100, g = 5, and various values of e30 between 0 and 0.001. It is the most frequent one. In this scenario, no capture into secondary resonances can be seen just before or after capture in the primary

FIG. 1. Tidal evolution in the present Mimas–Tethys i 1 i 3 resonance of the libration amplitude of the resonant angle (a), of the eccentricity of Tethys e3 in logarithmic scale (b), and of the inclinations of Mimas i 1 and Tethys i 3 (c), for n 01 = 2422.5072 rad year−1 , n 03 = 213.1732 rad year−1 , σ 0 = 180.64◦ , Q3 = 100, g = 5, and the following present values of e3 : e30 = 0 (curve 1), e30 = 0.000235 (curve 2), and e30 = 0.001 (curve 3). The dotted lines on Fig. 1b show the capture times into the i 1 i 3 resonance for the three curves. See also text and Fig. 2.

112

CHAMPENOIS AND VIENNE

from e30 = 0 (curve 1), e30 = 0.000235 (curve 2), and e30 = 0.001 (curve 3). The times of capture are tc ≈ −16.70 × 103 Qs years for trajectory 1, tc ≈ −15.79 × 103 Qs years for trajectory 2, and tc ≈ −14.90 × 103 Qs years for trajectory 3. The Mimas– Tethys resonance is thus found to be 11% younger with curve 3 than with curve 1. The results of our runs also reveal noticeable jumps in the mean value of the libration amplitude, at t ≈ −10.4 × 103 Qs years for e30 = 0.000235 and at t ≈ −10.2 × 103 Qs years for e30 = 0.001. The latter one is less important and is hidden by the oscillations of curve 3 on Fig. 1a. The variation of the eccentricity of Tethys corresponding to trajectories 1, 2, and 3 is plotted on Fig. 1b, as well as their corresponding time of capture in the i 1 i 3 resonance (dotted lines). A sharp downward variation may be seen on curve 1. Examining our results, we find that e3 = 8 × 10−5 on capture into the i 1 i 3 resonance according to trajectory 1, e3 = 0.002 according to trajectory 2, and e3 = 0.008 according to trajectory 3. We deduce from it that the condition e30 ∈ [0, 0.001] leads to

FIG. 2. Tidal evolution of the Mimas–Tethys system (the arrow indicates the direction of evolution) through the primary i 1 i 3 libration zone (V-shaped zone) and the main secondary resonances ( j/k curves; see text), for the three initial values e30 = 0 (1), e30 = 0.000235 (2), e30 = 0.001 (3), obtained from the runs shown on Fig. 1. Label “0” is for the present time. Inside the libration zone, r = 0.629308, and what is α13 oscillates around the exact resonance value α13 r . Outside shown on the figure are the borders of the libration of α13 around α13 the libration zone, the points represent real values of α13 .

e3 ∈ [8 × 10−5 , 0.008] at time of capture into the i 1 i 3 resonance. The variation of the inclinations of Mimas and Tethys since t = −18 × 103 Qs years is plotted in Fig. 1c. We notice that the value of i 1 before capture is all the higher as e30 is high. This is to be related to Fig. 1a. As for i 3 , it has been varying too little since capture time to be really sensitive to e3 . Hence only the value of i 1 before capture changes according to e30 : we have i 1 = 0.45◦ for e30 = 0 (curve 1), i 1 = 0.58◦ for e30 = 0.000235 (curve 2), and i 1 = 0.70◦ for e30 = 0.001◦ (curve 3). We find therefore the following interval of values for i 1 before capture: i 1 ∈ [0.45◦ , 0.70◦ ]. In order to get a better insight of the phenomena at work in the system, we have plotted on Fig. 2 the preceding numerical √ trajectories (from our complete equations) in an (α13 , γ1 γ3 ) plane. The arrow shows the direction of evolution, the “V” represents the libration zone, and the curves labeled j/k are the j/k secondary resonances, which arise when we have jω + k σ˙ = 0, where ω is the libration or circulation frequency of the i 1 i 3 resonance, regarded as a pendulum (the “V” and the curves for secondary resonances are both obtained by regarding the i 1 i 3 resonance as a pendulum; see CV for details). Note that outside the “V”, the points of the evolutionnary paths represent real values of α13 , while inside the “V”, the two paths represent only the borders of the libration of α13 around the exact resonance r value α13 . These paths are deduced from the results shown in Figs. 1a and 1c and from the following relation between α13 and φ˙ (see Appendix A in CV): r =− α13 − α13

1 2 3 φ˙ . 6 n3

(10)

Figure 2 shows that, when the three paths join, after capture in the primary resonance, they follow a trajectory close to the one that the system would follow if only the resonant argument of the primary i 1 i 3 resonance was considered in our complete equations (see Appendix). That is to say, they follow the “single resonance theory” also called “adiabatic invariant theory” (see CV for details). The jumps in the mean value of the libration amplitude (Fig. 1a) are related to the crossing of the 1/2 secondary resonance. In the same way, the downward variation seen on curve 1 of Fig. 1b is related to the crossing of the 1/1 secondary resonance. A second possible scenario has been obtained from the following parameters and initial conditions: Q3 = 100, g = 5, n 01 = 2422.5077 rad year−1 , n 03 = 1213.1732 rad year−1 , σ 0 = 180.64◦ , and e30 = 0.001. We have therefore changed very slightly the mean motion of Mimas compared to curve 3 of Fig. 1. The results, displayed on Figs. 3 and 4, reveal the possibility of a capture in the 1/2 secondary resonance between σ˙ and the libration frequency of the i 1 i 3 resonance: the capture in the primary resonance takes place sooner than in the first scenario (about −17.7 × 103 Qs years; Fig. 3a). Moreover, the value of i 1 before capture is lower (i 1 ≈ 0.18◦ ; Fig. 3b). As for the behavior of i 3 , it is unchanged compared to Fig. 1b. After capture, there are large oscillations in the libration amplitude (Fig. 3a), followed by a jump to about 30◦ . Then there is a smooth decrease in the amplitude (Fig. 3a), which follows the adiabatic invariant theory (Fig. 4) and the oscillations become larger as the 1/2 secondary resonance is approached: the libration amplitude falls down to about 5◦ (Fig. 3a). Once the system is trapped in the 1/2 secondary resonance (Fig. 4), the amplitude increases (from t ≈ −14.9 × 103 Qs years to t ≈ −10.3 × 103 Qs years; Fig. 3a). This phenomenon leads the system closer to the border of the i 1 i 3

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

113

FIG. 3. Tidal evolution in the present Mimas–Tethys i 1 i 3 resonance of the libration amplitude of the resonant angle (a), and of the inclinations of Mimas i 1 and Tethys i 3 (b), for n 01 = 2422.5077 rad year−1 , n 03 = 1213.1732 rad year−1 , σ 0 = 180.64◦ , Q3 = 100, g = 5, and e30 = 0.001. These figures reveal a capture into the 1/2 secondary resonance (see text and Fig. 4).

FIG. 5. Same as Fig. 1 for Q3 = 385, g = 3.85, n 01 = 2430.35 rad year−1 , n 03 = 1213.1743 rad year−1 , σ 0 = 0◦ , σ˙ = 2π/200.15 rad year−1 , and e30 = 0.001. These figures reveal a capture into the 1/1 secondary resonance (see text and Fig. 6).

FIG. 4. Same as Fig. 2, obtained from the run shown on Fig. 3. The system evolves until it gets trapped in the primary i 1 i 3 resonance (1) and then evolves in a chaotic way until it experiences a big variation that leads it very near the center of the island of this resonance (2). Then it evolves according to the adiabatic invariant theory (3), gets trapped in the 1/2 secondary resonance (4), quits it before the chaotic zone is reached (5), and evolves thereafter according to the single-resonance theory again, up to the present time.

resonance (Fig. 4). Afterward, the libration amplitude decreases again smoothly up to the present time (Fig. 3a). A third possible scenario, shown on Figs. 5 and 6, has been obtained with n 01 = 2430.35 rad year−1 , n 03 = 1213.1743 rad year−1 , σ 0 = 0◦ , Q3 = 385, g = 3.85, and e30 = 0.001, σ˙ being fixed throughout the integration to its current value (with period

114

CHAMPENOIS AND VIENNE

result obtained for e30 = 0.001. The method we used consists of restricting σ˙ to its secular and oblateness parts (see Eqs. (A.10) and (A.16)): dσ = dt

FIG. 6. Same as Fig. 2, obtained from the run shown on Fig. 5. After capture into the primary i 1 i 3 resonance, the system evolves until it gets trapped in the 1/1 secondary resonance (1). This leads it toward the chaotic zone of the primary resonance (2), in which the secondary resonance is disrupted. Afterward, it gets captured in the primary resonance again, across the 1/2 secondary resonance (3), and finally evolves up to the present time (4).

200.15 years: see the discussion below). The increase in the libration amplitude (Fig. 5a) is this time related to capture in the 1/1 secondary resonance, as revealed by Fig. 6 (the capture time on Figs. 5a–5c is not to be taken into account because of the value taken for n 01 ). Tidal evolution follows again the singleresonance theory before and after capture in that secondary resonance. The inclination of Mimas before capture is found to be quite low (about 0.03◦ ; Fig. 5c). As for the behavior of e3 , we notice that it grows abruptly during capture in the 1/1 secondary resonance, as can be seen on Fig. 5b. This figure thus displays a vanishingly low value for e3 on capture into the primary resonance. The peak in e3 is an interesting phenomenon which will be explained in Section 5. The main problem with this latter scenario lies in the fact that it has been obtained with a fixed σ˙ . In order to know if it is important to allow σ˙ to vary, we looked for the variation of σ˙ since capture time, for different values of e30 . Figure 7 shows the

FIG. 7. Tidal evolution of the frequency σ˙ for e30 = 0.001. See text for comments.

µ

dσ dt

¶

µ + S

dσ dt

¶ .

(11)

O

We notice on Fig. 7 that the period T of σ˙ has been slowly increasing since capture time up to the present value. Before capture, the increase was faster, due to the faster decrease in n 1 . The present value of T is found to be 204 years. There is therefore a slight discrepancy between the value thus obtained for the present time and the TASS value (200.15 years). Part of it is due to the method we used for determining σ˙ . As a matter of fact, TASS uses the frequency analysis method. When applied to our results, this methods yields of value closer to that of TASS (201.3 years). There remains a slight difference which comes from the truncation we had to do in our equations so as to speed up the computations: TASS equations are more complete than ours, insofar as they take account of the eight major saturnian satellites and the Sun and contain all critical terms (that is, secular, resonant, and solar terms) up to the ninth degree in eccentricity–inclination and second order in the masses and oblateness coefficients (third order for the J2 terms)). If we rely on the frequency analysis method, which does not neglect any term, Eq. (11) seems to overestimate T of about 3 years with respect to the frequency analysis method. As for the value of T on capture, our results show it is submitted to slight changes with e30 (due to the fact that the time of capture varies with e30 (top of this section)): according to Eq. (11), we have T = 187 years for e30 = 0.000235 and T = 190 years for e30 = 0.001. Finally we may say, relying on the frequency analysis method, that σ˙ has experienced a variation of about 8%, from about 185 years to the present value of 200 years. Hence we believe that the capture into the 1/1 secondary resonance obtained with a fixed σ˙ means that such a capture is still possible for a varying σ˙ (in fact, a shorter period tends to move the 1/1 secondary resonance location away from the separatrix, which decreases the risk it merges in the chaotic zone (see Section 5)). As a conclusion, our numerical integrations show that three main scenarii can be distinguished concerning the past evolution of the Mimas–Tethys system inside the main i 1 i 3 resonance: either the system is not trapped in any secondary resonance or it is trapped in the 1/1 or the 1/2 secondary resonance. In the first scenario, the i 1 i 3 resonance is encountered with i 1 ranging from 0.45◦ to 0.70◦ according to e3 , which ranges from 8 × 10−5 to 0.008. It should, however, be noted that the quoted values for e3 and i 3 depend strongly on the ratio Qs /Q3 : in our runs, we took it to be equal to 14,000/100 = 1400. If it should be higher (with for instance Qs = 14,000 and Q3 = 200), then e3 and i 1 on capture would both be weaker. They would be higher in the opposite case. In the second and third scenarii, the system passes some time in the 1/2 and 1/1 secondary resonances, which results in a lower value of i 1 on encounter (respectively 0.18◦ and 0.03◦ ).

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

However, the consequences on e3 of the trapping in the 1/2 and 1/1 resonances are not the same: concerning the second scenario, it has no effect at all: e3 behaves as if the system was not trapped in the 1/2 secondary resonance. This results in a rather high value of e3 on encounter of the i 1 i 3 resonance: e3 = 0.008. On the contrary, the trapping in the 1/1 resonance in the third scenario has a great effect on e3 , with the consequence that the value of e3 before capture is nearly zero. A second difference between the second and third scenarii concerns the value of the libration amplitude of the i 1 i 3 resonance on escaping from the secondary resonance: it is ≈110◦ for the second scenario, whereas it reaches ≈165◦ for the third one. The 1/1 secondary resonance thus leads the system far closer to the border of the i 1 i 3 resonance than the 1/2 secondary resonance. These capture scenarii move backward the time of capture into the major resonance, but their main interest is to call into question the result i 1 = 0.42◦ before capture found by Allan and to show that e3 may have been greater in the past. 5. SURFACES OF SECTION

˙ of the Mimas– Figure 8 shows the surfaces of section (2φ, 2φ) Tethys system for trajectories 2 (left) and 3 (right), shown in

115

Fig. 1, both at present time (up) and capture time (down), computed by considering the i 1 i 3 resonance as a pendulum perturbed by the external frequency σ˙ , considered constant (but we took two different values at capture time and present time; see CV for details). This figure reveals narrow chaotic borders of the i 1 i 3 resonance at present time, but large ones at capture time (when the eccentricity of Tethys is bigger), as well as islands associated with the 1/1 and 1/2 secondary resonances. In fact, other computations show that the chaotic zone is all the larger as e3 is high on capture. In the same way, the secondary resonance islands that can be seen on Figs. 8c and 8d (shown by arrows) become larger as e3 is increased. This explains why capture takes place later and later as e30 grows (Fig. 1). We may now also explain the differences that we noticed between the two jumps in the mean value of the libration amplitude, caused by the crossing of the 1/2 secondary resonance, mentioned in Section 4.2: the jump is bigger for e30 = 0.001 and shifted in time compared to the one associated with e30 = 0.000235, because the 1/2 secondary resonance island (in the ˙ is larger for e30 = 0.001 than for e30 = phase plane (2φ, 2φ)) 0.000235. Also we notice on curve 3 of Fig. 1a that there are large oscillations of the libration amplitude just after the crossing of the

FIG. 8. Surfaces of section of the Mimas–Tethys system at present time (top) and capture time (bottom) for e30 = 0.000235 (left) and e30 = 0.001 (right). The arrows indicate the secondary resonance islands. See text for comments.

116

CHAMPENOIS AND VIENNE

chaotic border of the i 1 i 3 resonance (the amplitude of the oscillations reach 18◦ ). As the system moves away from the chaotic border, the amplitude of the oscillations becomes smaller, especially after the crossing of the 1/2 secondary resonance. We see that the island associated with this resonance in the phase plane ˙ acts as a barrier to the chaotic diffusion. (2φ, 2φ) At last, the differences in the two scenarii involving a capture in a secondary resonance (scenarii 2 and 3) can now be better understood: in the third scenario, e3 increases sharply when the system is trapped in the 1/1 secondary resonance: this strengthens the island associated to this resonance in the phase plane, which becomes larger and larger. As this island becomes larger and larger, it moves toward the chaotic border of the i 1 i 3 resonance (Fig. 6). Therefore, the 1/1 secondary resonance island is not disrupted until it merges into the chaotic zone, which explains the high libration amplitude (≈165◦ ; Fig. 5a) on escape from this secondary resonance (Fig. 5a). This also explains the presence of a spike in e3 on Fig. 5b: shortly after escaping from the secondary resonance, the system escapes from the primary i 1 i 3 resonance and is captured again in it immediately after. Therefore it crosses the chaotic zone twice, in opposite directions, which results in a sharp increase in e3 , followed by a sharp decrease. On the contrary, in the second scenario, e3 is not increased as the system is trapped into the 1/2 secondary resonance. Hence the system moves toward the chaotic zone with an ever shrinking island, since e3 keeps decreasing because of tidal effects in Tethys (curve 3 on Fig. 1b). This entails a disruption of the secondary resonance before the chaotic zone is reached (whence the low amplitude on escape (≈120◦ ; Fig. 3a)). 6. CAPTURE PROBABILITY

6.1. Principle of Sinclair’s Method (1972) Sinclair (1972) considers Eqs. (6) and shows that capture into the i 1 i 3 resonance can be accounted for by the secular increase in the inclinations just before capture, which results in a decrease in the kinetic energy of the system and may therefore cause it to get captured if the circulation of the argument ϕ = 2φ was on the point of changing direction because of the tides. To evaluate the capture probability, he determines the value v of ϕ, ˙ as ϕ passes the unstable equilibrium point (near ϕ = −π), such that ϕ˙ = 0 when next at ϕ = π , without the system being captured between the two points. Then he determines the value u of ϕ˙ at ϕ = −π such that capture should just occur. Figure 9 shows on the (ϕ, ϕ) ˙ phase plane the trajectories such that ϕ˙ = u and ϕ˙ = v at ϕ = −π (it should be kept in mind that the axes ϕ = −π and ϕ = π are to be identified as the pendulum phase space is actually a cylinder). Assuming that all values of ϕ˙ 2 between 0 and v 2 are equally likely, which is equivalent to assuming that the system has random energy when at ϕ = −π, Sinclair then deduces the following formula for the probability of capture: u2 P= 2 = v

µ ¶2 u . v

˙ and trajectories of the Mimas–Tethys FIG. 9. Phase plane (ϕ = 2φ, ϕ˙ = 2φ) system such that ϕ˙ = u and ϕ˙ = v at ϕ = −π (the axes ϕ = −π and ϕ = π should be identified; see text).

To determine u and v, Sinclair resorts to numerical integrations. Integrating the same Eqs. (6) as Sinclair did, and using formula (12), we confirmed the value 0.04 he had found for the probability of capture of the Mimas–Tethys system into the i 1 i 3 resonance. 6.2. Generalization to Chaotic Systems 6.2.1. The method. We saw in Section 5 that the frequency σ˙ was responsible for a chaotic layer of noticeable width at the transition between libration and circulation of the i 1 i 3 resonance. In that case, the hyperbolic point no longer exists, and Sinclair’s method needs some modification. If we consider the system just before the circulation changes direction, there is now one more choice in addition to capture and inverse circulation: the system may also librate and circulate in turn in a chaotic way. The duration of this chaotic period depends on the initial condition chosen, but tidal effects place bounds on it, urging the system to choose for good between capture and inverse circulation. On inspecting formula (12), we notice that, in fact, it is equivalent to ¶2 µ l[0, u] , (13) P= l[0, u] + l[u, v] where l[0, u] is the length of the interval of values for which the system gets captured, and l[u, v] the length of the interval of values for which the circulation of the system changes direction. Accordingly, we denote, in case of chaotic systems, by m(C) the measure (on the set of real numbers) of the set C of values for which the system is captured for good and by m(E) the measure of the set E of values for which the circulation of the system changes direction (that is to say, escapes from resonance), we can thus define a probability by the following formula: µ

(12)

P=

m(C) m(C) + m(E)

¶2 .

(14)

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

A priori, the sets C and E are no more intervals. Now, two difficulties arise: the first one is to define the time length after which capture may be regarded as certain. Neishtadt (1997) chooses for a perturbation with amplitude ² the duration T = 1/². The duration corresponding to the weakest of our three perturbations being about 3000 years, we took T = 10,000 years out of security. The second problem lies in the way of numerically computing C and E. Our numerical experiments were made using the system we describe in the Appendix, analogous to the equations originally integrated by Sinclair, with in addition the major terms containing the angle σ . As C and E are no more intervals, we have to sweep the energy range at φ = π between the energy below which φ˙ is always negative and the energy above which φ˙ is always positive. In fact, we noticed that C and E can be written as [ C= Ci i

E=

[

(15)

Ei ,

i

where Ci and Ei are indeed intervals. So, by sweeping the energy range defined above, we determine the bounds of the major intervals Ci for which certain capture occurs, and of the major intervals Ei for which escape occurs. 6.2.2. Results. The probability found by Sinclair (0.04) seems to be weak. Of course, the system he integrated did not take account of the frequency σ˙ which Vienne and Duriez (1992) have detected since that time. Using the method described in Section 6.2.1, we can determine the probabilities corresponding to initial conditions on capture into the major resonance for various values of the eccentricity of Tethys on capture. We thus obtain a capture probability P that depends on e3 , given by Table III. P appears to vary strongly with e3 , since we have P = 0.033 for e3 = 0.0004 and P = 1 for e3 = 0.0016! We therefore see that the present uncertainty on e30 (that is to say, the uncertainty on e3 at capture time) results in a big uncertainty on the capture probability. However, it seems likely to be bigger than that determined by Sinclair. The method we developped also allows us to compute the overall capture probability of the two scenarii involving a capture TABLE III Probability of Capture into the Present i 1 i 3 Resonance as a Function of the Value e3 of Tethys’s Eccentricity on Capture e3

P

0 0.0004 0.0005 0.0008 0.0016

0.04 0.033 0.25 0.76 1

117

in the 1/1 and 1/2 secondary resonances, if only we can compute the probability of capture in these secondary resonances. This can be done in an analytical way for secondary resonances such as the 1/1 secondary resonance (see Malhotra 1990), but such calculations go beyond the scope of this paper. Furthermore, no model exists at the present time for the 1/2 secondary resonance. 7. CONCLUSION

We have investigated in this work the role of the 200-year-long period discovered by Vienne and Duriez (1992) on the dynamics of the Mimas–Tethys system presently trapped in the i 1 i 3 resonance at the 2 : 4 commensurability. This frequency σ˙ is in fact a ˙ 1− 3Ä ˙ +$ ˙ 3 ) whose influence on the secular frequency (σ˙ = 12 Ä 2 3 metric elements arises from three major arguments, whose terms in the expansion of the perturbative potential are all factored by the eccentricity e3 of Tethys. This latter one is badly known at present, but certainly lower than 0.001. We use simplified equations taking into account all the terms which generate this frequency (oblateness terms up to J23 and secular terms from the six inner classical saturnian satellites) and dissipation in Saturn and in Tethys (the latter one was not considered in previous studies). With these equations, we show from numerical integrations that, according to the value chosen for the present eccentricity e30 of Tethys, the presence of σ˙ entails two major differences in the evolutionary behaviour as described by Allan: first, the inclination i 1 of Mimas before capture in the primary resonance may be higher (if there is no capture in a secondary resonance) or lower (if there is a capture in the 1/1 or 1/2 secondary resonance) than that deduced by Allan (i 1 = 0.42◦ ). In the first case, assuming a ratio r = Qs /Q3 = 140 means i 1 ∈ [0.45◦ , 0.70◦ ] and e3 ∈ [8 × 10−5 , 0.008] before capture, according to e30 . In the second case, a capture in the 1/2 resonance (with r = 140) leads to i 1 ≈ 0.18◦ before capture, while a capture in the 1/1 resonance (with r = 36) gives i 1 ≈ 0.03◦ . Second, the eccentricity of Tethys may have been greater in the past, although a capture in the 1/1 secondary resonance would imply a vanishingly small value for e3 on capture in the primary resonance. Therefore, Allan’s scenario (1969) is called into question. However, it remains possible for very low values of e30 (a value e30 = 0 means i 1 = 0.45◦ before capture, which is near the value 0.42◦ found by Allan). Three main scenarii are therefore to be considered, according as there is no capture in a secondary resonance, capture in the 1/2 secondary resonance, or capture in the 1/1 secondary resonance. Hence calculation of the probability of capture is more complicated than given by Sinclair, all the more as we need to compute, in the second and third scenarii, probabilities of capture in secondary resonances. Furthermore, if e3 6= 0, the chaotic layer must be taken into account to get the probability of capture in the primary resonance. We developed a method of computing the probability in this last case, by generalizing Sinclair’s method (1972). At first, we assumed e3 to be identical to zero, and we confirmed Sinclair’s result P = 0.04. Then, considering

118

CHAMPENOIS AND VIENNE

e3 6≡ 0, we showed that P is extremely sensitive to e30 : it may be 0.03 as well as 1. Therefore, there exists an evolutionary path leading to the present state with reasonably high probability, characterized by relatively high values for e3 and i 1 at the time of encounter, but the possibility of several captures and escapes of the primary resonance (characterized with low values for i 1 at the time of encounter), due to the presence of secondary resonances, cannot be discarded. As a result, the time of capture in the i 1 i 3 resonance is uncertain, between t ≈ −17.7 × 103 Qs years and t ≈ −14.9 × 103 Qs years. At the end of this work, we realize that considering an eccentric orbit of Tethys upsets the vision of the dynamics of the Mimas–Tethys system we had so far. However, the law of the dependence of the probability of capture into the present i 1 i 3 resonance with the present eccentricity of Tethys still needs some investigation. APPENDIX: EQUATIONS OF MOTION Let R13 (resp. R31 ) be the disturbing function of Tethys on Mimas (resp. Mimas on Tethys). In the expressions R13 and R31 , we may content ourselves with the first order in the masses, since the masses of Mimas and Tethys are so small (see Table II), and set: 1 + m 1 ≈ 1 + m 3 ≈ 1. We thus have (see, e.g., Brouwer and Clemence 1960) µ

¶

r1 · r 3 1 a3 m3µ − a3 a3 113 |r3 |3 ¶ µ r1 · r 3 1 a3 , = m1µ − a3 a3 113 |r3 |3

R13 = R31

(A.1)

where µ = n 21 a13 = n 23 a33 = G MS , G is the gravitational constant, r1 and r3 are the coordinates of Mimas and Tethys, and 113 is the distance between them. In the neighborhood of the 2 : 4 i 1 i 3 resonance, the low-frequency terms are, up to the third degree in eccentricity–inclination, the terms with arguments φ, 3φ + σ, φ + σ , and −φ + σ . We recall that the arguments of the type Ä1 − Ä3 are short-period arguments when compared to 2φ and σ (their period is of the order of the year against more than 70 years for 2φ and about 200 years for σ ). The arguments φ, 3φ + σ , φ + σ , and −φ + σ are only present in the direct part of R13 and R31 (that is, a3 /113 ). So we have the following expressions for the L and R L of the disturbing functions R and R : long-period parts R13 13 31 31 ¶ µ 1 a3 m3µ a3 113 L ¶ µ 1 a3 = m1µ , a3 113 L

TABLE A.I Analytical Expressions of the Functions fi (α) for the Arguments 2φ, ψ + 2φ, ψ, ψ − 2φ, with Their Value for α0 = 0.63064 (TASS1.6) Argument

i

f i (α)

2φ

0

−αb3/2 (α)

ψ + 2φ

1

(4) d (4) 3αb3/2 (α) + 14 α 2 dα b3/2 (α)

ψ

2

d b3/2 (α) 2αb3/2 (α) + 12 α 2 dα

9.7082

3

(2) d (2) −αb3/2 (α) + 14 α 2 dα b3/2 (α)

0.22189

ψ − 2φ

a3 113

= f 0 (α13 )γ1 γ3 cos(2φ) + f 1 (α13 )γ12 e3 cos(3φ + σ ) L

+ f 2 (α13 )γ1 γ3 e3 cos(φ + σ ) + f 3 (α13 )γ32 e3 cos(−φ + σ ), (A.3) where ( f i )i=0,3 are functions of the classical Laplace’s coefficients and are given in Table A.I. In order to get the total gravitational disturbing function R1 that Enceladus, Tethys, Dione, Rhea, and Titan exert on Mimas, we have to add to L the secular parts (R S ) R13 1 j j=1,6 of the disturbing functions (R1 j ) j=1,6 of these

(2)

5.2379

(A.4)

We recall that these secular terms are only the terms whose arguments contain no angles, since the angle-dependent secular terms have been averaged out. In the L the secular parts (R S ) same way, we have to add to R31 3 j j=1,6 of the disturbing functions (R3 j ) j=1,6 of these satellites on Tethys to get the total gravitational perturbation on Tethys: L S S S S S R3 = R31 + R31 + R32 + R34 + R35 + R36 .

(A.5)

The actions of Japet, the Sun, and the small satellites of Saturn are not taken into account because of their weak effects in the generation of the perturbative frequency σ˙ (this was checked for Japet and the Sun with TASS1.6, which takes them into account). The effect of Saturn’s oblateness will be considered farther. Denoting by 1i j the distance between the satellites with numbers i and j and using the subscript S to denote the secular part, the expressions of R1S j and R3S j take the following form (see, e.g., Duriez 1989): ¶ µ aj 1 m jµ aj 11 j S ¶ µ 1 a3 = m jµ a3 13 j S ¶ µ aj 1 m jµ = aj 13 j S

R1S j = R3S j

if j = 1, 2

(A.6)

if j = 4, 5, 6.

When the secular part is restricted to third-degree in eccentricity–inclination and averaged on the angle-dependent terms, we have (with i < j)

(A.2)

¶

(2)

−1.6509

L S S S S S + R12 + R13 + R14 + R15 + R16 . R1 = R13

µ

with µ

(3)

satellites on Mimas, so that

L = R13 L R31

f i (α0 )

aj 1i j

¶

¢ ¢ ¡ ¡ = A(αi j ) + B(αi j ) ei2 + e2j + C(αi j ) γi2 + γ j2 ,

(A.7)

S

with A(αi j ) =

1 (0) b (αi j ) 2 1/2

B(αi j ) =

d (0) d 2 (0) 1 1 b1/2 (αi j ) + αi2j 2 b1/2 (αi j ) αi j 4 dαi j 8 dαi j

(A.8)

1 (1) C(αi j ) = − αi j b3/2 (αi j ). 2 The numerical values of the coefficients A, B, C are given in Table A.II for the values αi j which involve Mimas or Tethys. We then use Lagrange’s equations for the variables a1 , a3 , γ1 , γ3 , e3 , φ, and σ with the perturbative functions R1 and

119

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

TABLE A.II Values of αi j , A(αi j ), B(αi j ), C(αi j ), and d A(αi j )/dαi j for Every Pair (i, j) (i < j) Involving Mimas or Tethys i– j 1–2 1–3 1–4 1–5 1–6 2–3 3–4 3–5 3–6

αi j 0.78026 0.63064 0.49258 0.35283 0.15223 0.80824 0.78108 0.55948 0.24140

A(αi j ) 1.2473 1.1306 1.0706 1.0335 1.0059 1.2807 1.2482 1.0960 1.0151

B(αi j ) 1.3674 0.39952 0.15366 0.059940 0.0090810 1.8535 1.3790 0.23856 0.024460

C(αi j )

d A(αi j )/dαi j

−5.4695 −1.5581 −0.61463 −0.23976 −0.036324 −7.4138 −5.5159 −0.95425 −0.097840

1.0996 0.55451 0.33609 0.20479 0.078148 1.2978 1.1047 0.42538 0.12912

Note. The αi j are from TASS1.6.

µ

µ

dφ dt dσ dt

£ dn 1 = 3n 21 α13 m 3 2 f 0 (α13 )γ1 γ3 sin(2φ) + 3 f 1 (α13 )γ12 e3 sin(3φ + σ ) dt ¤ + f 2 (α13 )γ1 γ3 e3 sin(φ + σ ) − f 3 (α13 )γ32 e3 sin(−φ + σ )

µ

¶ = S

dλ1 dt

¶

µ µ ¶ µµ ¶ ¶ ¶ dλ3 1 dÄ3 dÄ1 −2 + + dt S 2 dt S dt S S (A.10)

µ µ ¶ µ ¶ ¶ 1 dÄ1 3 dÄ3 d$3 = − + 2 dt S 2 dt S dt S

¶ S

and µ

µ

µ

R3 , and we retain only the largest terms in the right-hand sides. The variables that we consider are free from short-period effects. We should therefore add to Lagrange’s equations second-order terms in the masses (see, e.g., Vienne (1991)) but we neglect them on account of their weakness (the mass of the most massive satellite, Titan, is less than M S /4000 (Table II)), and we use the same notation for both types of variables. We obtain:

µ

µ

dλ1 dt dλ3 dt

dÄ1 dt dÄ3 dt d$3 dt

¶ = −2n 1 S

"

¶ = 2n 3 S

¶ S

i=2

d A(α1i ) dt

# ¶ X 6 d A(αi3 ) 2 d A(α3i ) αi3 + A(αi3 ) + m i α3i dαi3 dα3i i=4

2 µ X i=1

=

n3 2

"

2 X

"

¶ = 2n 3 S

2 m i α1i

6 n1 X = m i α1i C(α1i ) 2 i=2

¶ S

6 X

(A.11)

m i C(αi3 ) +

i=1 2 X

6 X

# m i α3i C(α3i )

i=4

m i C(αi3 ) +

i=1

6 X

# m i α3i B(α3i ) .

i=4

The values of d A(αi j )/dαi j involving Mimas or Tethys are given in Table A.II for the corresponding values αi j given by TASS1.6. The variation of (dλ1 /dt)S , (dλ3 /dt)S , (dÄ1 /dt)S , (dÄ3 /dt)S , and (d$3 /dt)S with n 1 and n 3 during capture in resonance (less than 1% (CV)) leads to negligible changes in n 1 , n 3 , (dÄ1 /dt), (dÄ3 /dt), and (d$3 /dt), so they are neglected. We then take for the secular variations of λ1 , λ3 , Ä1 , Ä3 and $3 the following constants (in rad year−1 ), issued from Tables I, II, and A.II

£ dn 3 = −6n 23 m 1 2 f 0 (α13 )γ1 γ3 sin(2φ) + 3 f 1 (α13 )γ12 e3 sin(3φ + σ ) dt ¤ + f 2 (α13 )γ1 γ3 e3 sin(φ + σ ) − f 3 (α13 )γ3 2 e3 sin(−φ + σ ) dγ1 1 = n 1 α13 m 3 [γ3 f 0 (α13 ) sin(2φ) + 2γ1 e3 f 1 (α13 ) sin(3φ + σ ) dt 4 + γ3 e3 f 2 (α13 ) sin(φ + σ )]

µ µ

1 dγ3 = n 3 m 1 [γ1 f 0 (α13 ) sin(2φ) − γ1 e3 f 2 (α13 ) sin(φ + σ ) dt 4 − 2γ3 e3 f 3 (α13 ) sin(−φ + σ )] £ de3 = n 3 m 1 γ12 f 0 (α13 ) sin(3φ + σ ) + γ1 γ3 f 2 (α13 ) sin(φ + σ ) dt ¤ + γ32 f 3 (α13 ) sin(−φ + σ )

with

µ µ (A.9) µ

dλ1 dt dλ3 dt

dÄ1 dt dÄ3 dt

¶ = −0.0050046 ¶

S

= −0.0078503 S

¶

= −0.0047843

(A.12)

S

¶

= −0.010654 S

¶

¶ ·µ dφ γ3 γ1 1 f 0 (α13 ) cos(2φ) + n3 m 1 = n 1 − 2n 3 + n 1 α13 m 3 dt 8 γ1 γ3 + 2n 1 α13 m 3 e3 f 1 (α13 ) cos(3φ + σ ) ¶ µ γ3 γ1 e3 f 2 (α13 ) cos(φ + σ ) + n 1 α13 m 3 − n3 m 1 γ1 γ3 ¸ µ ¶ dφ + 2n 3 m 1 e3 f 3 (α13 ) cos(−φ + σ ) + dt S

We notice that Eqs. (A.9) are singular in φ and σ , so we change to the following variables, in order to have regular equations:

¸ · dσ γ3 γ1 1 f 0 (α13 ) cos(2φ) − 3n 3 m 1 = n 1 α13 m 3 dt 8 γ1 γ3 ¸ · m 1 n 3 γ12 n 1 α13 m 3 e3 f 1 (α13 ) cos(3φ + σ ) + + 4 e3 · · ¸¸ n 1 α13 m 3 γ3 e3 3e3 γ3 + + n 3 m 1 γ1 − f 2 (α13 ) cos(φ + σ ) 8γ1 e3 8γ3 ¸ µ ¶ · m 1 n 3 γ32 dσ −3m 1 n 3 e3 f 3 (α13 ) cos(−φ + σ ) + , + + 4 e3 dt S

We now have four degrees of freedom instead of 3 2 , but this is not really a problem, since we intend to do numerical integrations. Let us point out that the consideration of w and u instead of σ and φ settles the numerical problems which may arise with Eqs. (A.9) when e3 is small or equal to zero, and when i 1 or i 3 is small. However, problems will subsist with the regular equations if i 1 or i 3 is strictly equal to zero, because the computer will have to deal with divisions of the type “zero divided by zero.” However, this is really a particular case, and, actually, variables (13) are adequate for our purposes. We must now add to Eqs. (A.9) (made regular) the terms due to the oblateness of Saturn. Saturn’s gravitational momenta are quite important (see Table II), so that we have, in order to get the full variations of the mean longitudes, nodes

d$3 dt

n 1 , n 3 , γ1 , γ3 , w = γ1 γ3 e

= 0.010594. S

√ −1φ

, w, ¯ u = γ1 γ3 e3 e 1

√ −1 σ

¯ , u.

(A.13)

120

CHAMPENOIS AND VIENNE

and pericenters due to the oblateness, to take account of the lowest-degree terms with J2 , J4 , J6 , J22 , J2 J4 , and J23 as a factor (see, e.g., Vienne (1991) for a description of the method which allows to get these terms; we simply note here that the J22 , J2 J4 , and J23 terms arise from the fact that we consider variables free from short-period oscillations). In the case of J2 , we also have to take into account the next degree in eccentricity–inclination. For the other small factors, they are negligible. Only for the J2 , J4 , and J22 terms is the variation with n (and with the eccentricities and inclinations for J2 ) nonnegligible. The other terms are thus taken constant. We then get µ

µ

µ

dÄ1 dt

dÄ3 dt

d$3 dt

µ

µ

dλ1 dt

dλ3 dt

¶ O

¶ O

¡ ¢ 3 7/3 = − J2 n 1 µ−2/3 ae2 1 + 2e12 − 2γ12 2 µ ¶ 15 45 2 11/3 + n 1 µ−4/3 J4 − J2 − 0.0037064 4 8

= O

¡ ¢ 3 7/3 J2 n 3 µ−2/3 ae2 1 + 2e32 − 2γ32 2 µ ¶ 15 63 2 11/3 + n 3 µ−4/3 − J4 + J2 + 0.000032293 4 8

O

µ ¶ 7 7/3 = 3J2 n 1 µ−2/3 ae2 1 + e12 − 7γ12 4 µ ¶ 15 160 2 11/3 + n 1 µ−4/3 − J4 + J2 + 0.0048096 4 11

O

µ ¶ 7 7/3 = 3J2 n 3 µ−2/3 ae2 1 + e32 − 7γ32 4 µ ¶ 15 160 2 11/3 J2 + 0.00015183, + n 3 µ−4/3 − J4 + 4 11

¶

¶

Allan, R. R. 1969. Evolution of Mimas–Tethys commensurability. Astron. J. 74, 497–506. Brouwer, D., and G. M. Clemence 1960. Methods of Celestial Mechanics. Academic Press, New York. Campbell, J. K., and J. D. Anderson 1989. Gravity field of the saturnian system from Pioneer and Voyager tracking data. Astron. J. 97, 1485–1495.

(A.14)

µ

(A.15)

with µ

µ ¶ µµ ¶ ¶ ¶ dλ3 1 dÄ3 dÄ1 + + dt O 2 dt O dt O O O Ã µ µ ¶ µ ¶ ¶ ¶ dσ 1 dÄ1 3 dÄ3 d$3 = − + . dt O 2 dt O 2 dt O dt O dφ dt

µ

¶

=

dλ1 dt

µ

¶

Champenois, S., and A. Vienne 1999. Chaos and secondary resonances in the Mimas–Tethys system. Celest. Mech. Dynam. Astron., in press. Dermot, S. F., R. Malhotra, and C. D. Murray 1988. Dynamics of the uranian and saturnian systems: A chaotic route to melting miranda? Icarus 76, 295–334. Dourneau, G. 1987. Observations et e´ tude du mouvement des huit premiers satellites de Saturne, Th`ese, Bordeaux. Duriez, L. 1989. Le developpement de la fonction perturbatrice. In Les Methodes Modernes de la Mecanique Celeste (D. Benest and C. Froeschle, Eds.), pp. 35– 62. Goutelas, France. Duriez, L., and A. Vienne 1997. Theory of motion and ephemerides of Hyperion. Astron. Astrophys. 324, 366–380. Gavrilov, S. V., and V. N. Zharkov 1977. Love numbers of the giant planets. Icarus 32, 443–449. Harper, D., and D. B. Taylor 1993. The orbits of the major satellites of Saturn. Astron. Astrophys. 268, 326–349.

where the subscript “O” is for “oblateness.” The literal expressions are from (Levallois and Kovalevski 1971). In these expressions, µ = G MS , J2 and ae are obtained from Table II, and e1 is taken from Table I. The numerical expressions are issued from the program developed by L. Duriez for the expansion of the right-hand side of Lagrange’s equations: each constant is the sum of the three contributions of J6 , J2 J4 , and J23 , expressed in rad year−1 . The oblateness terms we must add to (regularized) Eqs. (A.9) are therefore: µ ¶ ¶ √ dw dφ = −1w dt O dt O µ ¶ µ ¶ √ du dσ = −1u , dt O dt O

This research was carried out thanks to the financial contribution of the PNP (Planetology National Program) of the National Institute of the Sciences of the Universe (INSU). We are very grateful to J. Henrard and A. Lemaˆıtre for fruitful discussions. We also thank L. Duriez for his careful reading of the manuscript and the referees A. T. Sinclair for his suggestions and S. Ferraz-Mello and T. Yokoyama for a very detailed report.

REFERENCES

¡ ¢ 3 7/3 = − J2 n 3 µ−2/3 ae2 1 + 2e32 − 2γ32 2 µ ¶ 15 45 2 11/3 J4 − J2 − 0.00011694 + n 3 µ−4/3 4 8

¶

ACKNOWLEDGMENTS

−2

(A.16)

Regularized Eqs. (A.9), to which we add the oblateness terms (A.15) and tidal dissipation on the mean motions of Mimas and Tethys and on the eccentricity of Tethys (see Section 3), form the system of equations we used for the numerical integrations presented in this paper.

Henrard, J. 1982. The adiabatic invariant: Its use in celestial mechanics. In Applications of Modern Dynamics to Celestial Mechanics (V. Szebehely, Ed.), pp. 153–171. Reidel, Dordrecht. Henrard, J., and A. de Vleeschauwer 1988. Sweeping through a second order resonance. Celest. Mech. Dynam. Astron. 43, 99–112. Kaula, W. M. 1964. Tidal dissipation by solid friction and the resulting orbital evolution. Rev. Geophys. Space Phys. 2, 467–539. Kozai, Y. 1957. On the astronomical constants of saturnian satellites system. Ann. Tokyo Astron. Obs. (Second Series) V-2, 73–106. Laskar, J. 1990. The chaotic motion of the Solar System: A numerical estimate of the size of the chaotic zones. Icarus 88, 266–291. Levallois, J.-J., and J. Kovalevsky 1971. G´eod´esie g´en´erale: G´eod´esie spatiale, Tome IV. Editions Eyrolles, Paris, France. Malhotra, R. 1990. Capture probabilities for secondary resonances. Icarus 87, 249–264. Malhotra, R. 1991. Tidal origin of the Laplace resonance and the resurfacing of Ganymede. Icarus 94, 399–412. Malhotra, R., and S. F. Dermott 1990. The role of secondary resonances in the orbital history of Miranda. Icarus 85, 440–480. Nakamura, T., and O. Abe 1977. Internal friction of snow and ice at low frequency. In Proceedings of 6th International Conference on Internal friction and Ultrasonic Attenuation in Solids (R. R. Hasiguti and M. Mikoshiba, Eds.), pp. 285–289. Univ. of Tokyo Press, Tokyo. Neishtadt, A. 1997. Scattering by resonances. Celest. Mech. Dynam. Astron. 65, 1–20. Peale, S. J., P. Cassen, and R. T. Reynold 1980. Tidal dissipation, orbital evolution, and the nature of Saturn’s inner satellites. Icarus 43, 65–72.

SECONDARY RESONANCES IN THE MIMAS–TETHYS SYSTEM

121

Showman, A., and R. Malhotra 1997. Tidal evolution into the Laplace resonance and the resurfacing of Ganymede. Icarus 127, 93–111.

Vienne, A., and L. Duriez 1992. A general theory of motion for the eight major satellites of Saturn. Astron. Astrophys. 257, 331–352.

Sinclair, A. T. 1972. On the origin of the commensurabilities amongst the satellites of Saturn. Mon. Not. R. Astron. Soc. 160, 169–187.

Vienne, A., and L. Duriez 1995. TASS1.6: Ephemerides of the major saturnian satellites. Astron. Astrophys. 297, 588–605.

Tittemore, W., and J. W. Wisdom 1988. Tidal evolution of the uranian satellites. I. Passage of Ariel and Umbriel through the 5:3 mean motion commensurability. Icarus 74, 172–230.

Vienne, A., L. Duriez, and S. Champenois 1996. The Mimas–Tethys and Enceladus–Dione systems. In Dynamics, Ephemerides and Astrometry in the Solar System (S. Ferraz-Mello, B. Morando, and J. E. Arlot, Eds.), pp. 143– 144. IAU Symposium 172, Paris, July 1995.

Tittemore, W., and J. W. Wisdom 1989. Tidal evolution of the uranian satellites. II. An explanation for the anomalously high orbital inclination of Miranda. Icarus 78, 63–89. Vienne, A. 1991. Th´eorie Analytique des Satellites de Saturne. Th`ese, Universit´e des Sciences et Techniques de Lille-Flandres-Artois.

Vienne, A., J. M. Sarlat, and L. Duriez 1992. About the secular acceleration of Mimas. In Chaos, Resonance and Collective Dynamical Phenomena in the Solar System (Ferraz-Mello, Ed.), pp. 219–222. IAU Symposium 152, Angra Dos Reis, Brazil, July 1991.