Frequency Map and Global Dynamics in the Solar System I

Icarus 152, 4–28 (2001) doi:10.1006/icar.2000.6576, available online at http://www.idealibrary.com on

Frequency Map and Global Dynamics in the Solar System I Short Period Dynamics of Massless Particles P. Robutel and J. Laskar Astronomie et Syst`emes Dynamiques, Institut de M´ecanique C´eleste, CNRS UMR 8028, 77 Av. Denfert-Rochereau, 75014 Paris, France E-mail: [email protected] Received February 28, 2000; revised November 15, 2000

Wisdom 1993, Levison and Duncan 1993) as this now becomes possible if one uses a symplectic integrator and large step size (Wisdom and Holman 1991, Levison and Duncan 1994). On the other hand, the dynamics of asteroids has now been investigated in great detail (Wisdom 1983; Milani and Knezevic 1990; Moons and Morbidelli 1993, 1995; Ferraz-Mello 1994, Nesvorny and Morbidelli 1998), and the analytical and semianalytical methods developed within the framework of the asteroid belt have been applied to study the Kuiper belt dynamics (Morbidelli et al. 1995). Even particles in the inner Solar System have been numerically integrated in order to search for possible stable zones (Evans and Tabachnik 1999). We have thus several fragmented visions of the dynamics of massless particles in the Solar System, which is understandable, as for each specific cases, numerical integrations are conducted over extremely long times. In the first part of this paper, we present a unified vision of the dynamics of particles in the Solar System, obtained by the construction of a complete cartography of the Solar System dynamics, using Laskar’s Frequency map analysis (FMA) (Laskar, 1990, 1993; Laskar et al. 1992). The main advantages of this approach, beside having a rigorous support (Laskar 1999), is to permit very short integration times, and to identify easily the various resonances involved. This approach will allow us to limit the integration times to a few million of years, and we will thus be able to include a huge amount of initial conditions, which will provide a complete view of single-particle dynamics in the Solar System. In a first stage, which is presented here, we will limit ourselves to the consideration of short-term dynamics, or more precisely to diffusion and chaotic behavior resulting essentially from short-period resonances. We will also provide, in Section 4, details on the geometry of the involved resonances, which has not been investigated very much so far (Morbidelli et al. 1995, Malhotra 1996; Morbidelli 1997). This will be applied in Section 6 to the determination, for 62 Kuiper belt objects, of a picture of the dynamics of their vicinity, which in some sense provides for these objects a dynamical ID, which can be used, even if their orbital elements are not yet perfectly known.

Frequency map analysis (FMA) is a refined numerical method based on Fourier techniques that provides a clear representation of the global dynamics of multidimensional systems, which is very effective for systems of 3 degrees of freedom and more, and was applied to a large class of dynamical systems (Solar System, galaxies, particle accelerators, etc.). FMA requires only a very short integration time to obtain a measure of the diffusion of the trajectories, and makes it possible to identify easily the location of the main resonances. We have performed a complete analysis of massless particles in the Solar System, from Mercury (0.38 AU) to the outer parts of the Kuiper belt (90 AU), for all values of the eccentricities, and several values for the inclinations. This provides a complete dynamical map of the Solar System, which is, in this first step, restricted to mean motion resonances. The precise extent of the resonant islands in the phase space has been determined, and the vicinity of the phase space has been plotted for the 62 best known Kuiper belt objects. We have thus set up a system for the analysis of the numerous planetary systems that are expected to be discovered in the near future. As an example, we present the application of this method to the understanding of the dynamics of the newly discovered v-Andromedae system. °c 2001 Academic Press Key Words: Solar System dynamics; resonances; Kuiper belt objects; asteroids dynamics.

1. INTRODUCTION

In the past decade, the understanding of the dynamics of celestial bodies in the Solar System has increased a great amount. This is largely due to the improvement of computer speed, which now makes it possible to track the evolution of planets or massless particles over durations of several hundred of million of years. Several groups have been conducting numerical experiments on the dynamics of massless particle in the outer Solar System and beyond, starting with Torbett (1989) and Gladman and Duncan (1990). These studies are motivated by the search for the understanding of the dynamics and long term evolution of the Kuiper belt objects (Luu 1994). For the recent studies, numerical integrations of massless particles in the outer Solar System are conducted over several billion of years (Holman and 4 0019-1035/01 $35.00 c 2001 by Academic Press Copyright ° All rights of reproduction in any form reserved.

5

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

More generally, we are setting up tools that make it possible to investigate rapidly the stability and global dynamics of any planetary system. As an application, we present in the last part an analysis of the dynamics in the newly discovered v-Andromedae system. 2. FREQUENCY MAP ANALYSIS

In this section, we provide a brief description of the FMA method. More details can be found in Laskar (1999) and references therein. 2.1. Frequency Maps According to the KAM theorem (see Arnold (1988)), in the phase space of a sufficiently close to integrable conservative system, many invariant tori will persist. Trajectories starting on one of these tori remain on it thereafter, executing quasiperiodic motion with a fixed frequency vector dependent only on the torus. The family of tori is parameterized over a Cantor set of frequency vectors, while in the gaps of the Cantor set chaotic behavior can occur. Although the frequencies are strictly speaking only defined and fixed on these tori, the frequency analysis algorithm will numerically compute over a finite time span a frequency vector for any initial condition. On the KAM tori, this frequency vector will be a very accurate approximation of the actual frequencies, while in the chaotic regions, it will provide a natural interpolation between these fixed frequencies. Let us consider an n-DOF Hamiltonian system close to integrable in the form H (I, θ) = H0 (I ) + ε H1 (I, θ), where H is real analytic for (I, θ) ∈ B n × Γn , B n is a domain of Rn , and Γn is the n-dimensional torus. For ε = 0, the Hamiltonian reduces to H0 (J ) and is integrable. The equations of motion are then for all j = 1, . . . , n ∂ H0 (I ) = ν j (I ), I˙j = 0, θ˙ j = ∂Ij

(1)

of the form |(k, ν)| > κε /|k|m

(3)

for which the perturbed system still possesses smooth invariant tori with linear flow (the KAM tori). Moreover, according to P¨oschel (1982), there exists a diffeomorphism 9 : Γn × Ä → Γn × B n , (ϕ, ν) → (θ, I ),

(4)

which is analytical with respect to ϕ and C ∞ in ν and on Γn × Äε transforms the Hamiltonian equations into ν˙ j = 0, ϕ˙ j = ν j . For frequency vectors (ν) in Äε , the solution lies on a torus and is given in complex form by its Fourier series z j (t) = z j0 eiν j t +

X

am (ν)eit ,

(5)

m

where the coefficients am (ν) depend smoothly on the frequencies (ν). If we fix θ ∈ Γn to some value θ = θ0 , we obtain a frequency map on B n defined as Fθ0 : B n → Ä, I → p2 (9 −1 (θ0 , I )),

(6)

where p2 is the projection on Ä ( p2 (φ, ν) = ν). For sufficiently small ε, the torsion condition (2) ensures that the frequency map Fθ0 is a smooth diffeomorphism. 2.2. Quasiperiodic Approximations FMA relies heavily on the observation that when a quasiperiodic function f (t) in the complex domain C is given numerically, it is possible to recover a quasiperiodic approximation of f (t) in a very precise way over a finite time span [−T, T ], several orders of magnitude more precisely than by simple Fourier analysis. Indeed, let X

f (t) = eiν1 t +

ak eihk,νit , ak ∈ C

(7)

k∈Z −(1,0...,0) n

which gives z j (t) = z j0 eiν j t in the complex variables z j = I j exp iθ j , where z j0 = z j (0). The motion in phase space takes place on tori, products of true circles with radii I j = |z j (0)|, which are described at constant velocity ν j (I ). If the system is nondegenerate, that is, if µ

∂ν(I ) det ∂I

¶

µ

∂ 2 H0 (I ) = det ∂I2

¶ 6= 0,

(2)

the frequency map F: B n → Rn ; (I ) → (v) is a diffeomorphism on its image Ä, and the tori are as well described by the action variables (I ) ∈ B n or in an equivalent manner by the frequency vector (ν) ∈ Ä. For a nondegenerate system, the KAM theorem still asserts that for sufficiently small values of ε, there exists a Cantor set Ä of values of (ν), satisfying a Diophantine condition

be a KAM quasiperiodic solution of an Hamiltonian system in B n × Γn , where the frequency vector (ν) satisfies a Diophantine condition (3). The frequency analysis NAFF will P N algorithm 0 ak0 eiωk t of f (t) from its provide an approximation f 0 (t) = k=1 numerical knowledge over a finite time span [−T, T ]. The frequencies ωk0 and complex amplitudes ak0 are computed through an iterative scheme. In order to determine the first frequency ω10 , one searches for the maximum amplitude of φ(σ ) = h f (t), eiσ t i where the scalar product h f (t), g(t)i is defined by h f (t), g(t)i =

1 2T

Z

T

f (t)g¯ (t)χ(t) dt,

(8)

−T

and where χ(t) is a weight function, that is, a positive and even RT function with 1/2T −T χ(t) dt = 1 (in our computations, we

6

ROBUTEL AND LASKAR

used the Hanning window filter, that is, χ1 (t) = 1 + cos(π t/T )). 0 Once the first periodic term eiω1 t is found, its complex am0 plitude a1 is obtained by orthogonal projection, and the process is restarted on the remaining part of the function f 1 (t) = 0 f (t) − a10 eiω1 t . It is also necessary to orthogonalize the set of 0 iωk0 t functions (e )k , when projecting f iteratively on these eiωk t . For a KAM solution, this algorithm allows a very accurate determination of the frequencies over the time span [−T, T ], which converges toward their true values as T increases. Indeed, for a weight function of the form

properties for a large class of dynamical systems: Solar System (Laskar 1990; Laskar and Robutel 1993), particle accelerators (Dumas and Laskar 1993; Laskar and Robin 1996), galactic dynamics (Papaphilippou and Laskar 1996, 1998) and standard Maps (Laskar et al. 1992; Laskar 1993). 2.4. FMA of a Massles Particle

The system formed by the Sun, N planets, and a massless particle is (after the reduction of the center of mass) a system with 3N + 3 degrees of freedom, and its study is difficult, as it involves a large number of parameters. The phase space is of di2 p ( p!)2 p (1 + cos π t) (9) mension 6N + 6, and the frequency space of dimension 3N + 3, χ p (t) = (2 p)! which is still very large for the Solar System with N = 8 (Pluto is excluded as it is considered as a massless particle), but if we have the following (Laskar 1999): we are interested only in the behavior of the massless partiPROPOSITION. For a KAM solution f (t) of the form (7), and cle, the problem becomes much simpler, as this particle does using the weight function χ (t) = χ p (t), the application of the not perturb the planetary motions. The frequencies correspondfrequency analysis algorithm over the time span [−T, T ], as ing to the planets will thus be fixed and it is sufficient to study described above, provides a determination ν1T of the frequency the frequency map (11) in a 3-dimensional space of action-like ν1 , which converges towards ν1 with the asymptotic expression variables. For any particle, the three action-like variables are the usual elliptical elements (a, e, i) (semi-major axis, eccentricity, for T → +∞ and inclination), while the three angle variables are the associµ ¶ p+1 2 p 2X ated angles (M, ω, Ä) (mean anomaly, argument of perihelion, (−1) 1 π ( p!) <(a ) k cos(Äk T ) + o 2 p+2 ν1 − ν1T = 2 p+1 2 p+2 and longitude of the node). ApT T k Äk In order to construct the frequency map, we will thus fix for (10) t = 0 all the initial angles (M0 , ω0 , Ä0 ) of the particle, and construct numerically the frequency map from the initial action Pp 2 with Äk = hk, νi − ν1 , A p = − π22 ( π6 − k=1 k12 ), and where variable (a0 , e0 , i 0 ) to the numerically determined fundamental <(x) is the real part of x. In particular, the use of a Hanning frequencies of the trajectories, which are the frequencies redata window ( p = 1) ensures that for a KAM solution, the ac- lated to mean motion n, perihelion frequency g, and node frecuracy regarding the determination of the main frequencies will quency s: be proportional to 1/T 4 , instead of 1/T for an ordinary FFT. Thus, the frequency analysis will easily allow the recovery of (12) F T : (a0 , e0 , i 0 ) → (n, g, s). the frequency vector (ν1 , ν2 , . . . , νn ). 2.3. FMA With the previous algorithm, it is possible to construct numerically a frequency map in the following way: (a) Fix all values of the angles θ j = θ j0 . (b) For any value J0 of the action variables, integrate numerically the trajectories of initial condition (J0 , θ0 ) over the time span T . (c) Search for a quasiperiodic approximation of the trajectory with the algorithm of the previous section. (d) Identify the fundamental frequencies ν of this quasiperiodic approximation. The frequency map is then FθT0 : J 7→ ν

(11)

and from the previous section, we know that for T → +∞, on the set of regular KAM solutions, FθT0 → Fθ0 . In particular, FθT0 should converge toward a smooth function of the set of invariant KAM tori. The destruction of the invariant tori, and thus the chaotic regions, can then be identified by the nonregularity of the frequency map FθT0 . Thanks to its precision this method was used to initiate studies of global stability and refined diffusion

In fact, in the present study, we will even restrict ourselves to the determination of the short-period frequency (mean motion), which will allow us to use very short integration times. Moreover, we will fix a specified value of the initial inclination i 0 , and will thus consider only the map FiT0 : (a0 , e0 ) → n.

(13)

It should be noted that even if the plane (a0 , e0 ) is not of codimension one (it does not cross all the trajectories of the phase space) it is dynamically representative in the sense that it crosses all resonances. Therefore the choice of other initial phases will not qualitatively affect our map (see Section 4): only the local width of the encountered resonances will be affected but not their dynamical effect. The trajectories of massless bodies are computed using the symplectic integrator of second-order MVS of the Swift package (Wisdom and Holman 1991, Levison and Duncan 1994), with a time step of 0.7 years for the inner Solar System, and forty times less for the inner system. These integrators make it

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

possible to simultaneously integrate several massless bodies. The integration time 2T is quite short as it is typically 2 Ma (see Section 3). Moreover, particles are removed from the computation under special conditions. During the integration, two criteria are checked at each time step. We stop the integration if the particle follows an ejection path (the eccentricity of the instantaneous Keplerian heliocentric orbit becomes greater or equal to 1), or if its semi-major axis becomes larger than 1000 AU. We stop also the computation in case of a close encounter with a planet, what corresponds to a distance inferior to a third of the Hill radius of the considered planet. Although the test particle could survive a close approach without ejection, or injection in the inner Solar System (Levison and Duncan 1997, Morbidelli 1997), its elliptic elements would be drastically altered and the trajectory would be extremely far from quasiperiodic; it is thus considered as strongly irregular. 2.5. Diffusion For particles that survive the integration, the quantity z(t) = a(t) exp(iλ(t)) is analyzed by the frequency analysis algorithm, which provides a quasiperiodic approximation over the time span [0, T ]

z 0 (t) = α0 exp(iν0 t) +

N X

α j exp(iν j t).

(14)

j

If the motion of the particle were Keplerian, the decomposition (14) would contain only one term: α0 exp(iν0 t) where |α0 | = a(0) and ν0 is the particle mean motion. As it is not the case, the quasiperiodic decomposition (14) contains several periodic terms. The amplitude |α0 | of the first term in (14) is close to the initial semi major axis of the particle, and the frequency ν0 to the initial mean motion, while the amplitudes of the other terms are from 100 to 1000 times smaller (except for strongly irregular trajectories). If the trajectory is quasiperiodic (which implies its stability for all times) the frequency ν0 is an invariant of the motion (this is not the case for the osculating mean motion). This frequency, which we will call “proper mean motion” and denote by n, is one of the fundamental frequencies that parametrize the KAM torus on which lies the trajectory. The two other frequencies of the trajectory are the “proper” precession frequencies of perihelion and node (g, s), but at this stage, we will not compute them, which would require much longer integration times. In order to measure the lack of regularity of an orbit, we study the diffusion in frequency space. The proper mean motion is computed first on the time interval [0, T ], which gives the frequency n (1) , while the quantity n (2) is obtained on [T, 2T ]. On a quasiperiodic solution, the two frequencies are equal, and if it is not the case, the quantity σ = 1 − n (2) /n (1) gives an estimation of the diffusion rate over the time T (Laskar 1993). This frequency diffusion rate can be translated into semi-major axis

7

diffusion using the Kepler law n 2 a 3 = µ as 1n 3 1a ≈− . n 2 a

(15)

2.6. Mean Motion Resonances Apart from the estimation of the diffusion, the determination of the proper mean motion allows us to locate very accurately the libration areas associated with the mean motion resonances. All the massless bodies that belong to the same planetary orbital resonance have the same proper mean motion n. In a ( p1 : p2 ) resonance with the jth planet, n must satisfy the relation p1 n j + p2 n = 0 where n j is the proper mean motion of the considered planet. This property is illustrated in Fig. 1, which shows the proper mean motion of 200 particles with initial semi-major axis a0 ranging from 40 to 50 AU, and orbiting around the Sun in the perturbation field of the four giant planets of the Solar System. The initial elliptic elements of the test particles are chosen on a straight line of the phase space (i 0 = M0 = ω0 = Ä0 ; e0 = 0.15). The line of initial conditions crosses many islands of resonance with Neptune. These resonances are characterized by the flatness of the frequency curve (Laskar et al. 1992, Laskar 1993), while outside of these regions, the slope of the curve is √ close to a Keplerian one, that is, dn/da ≈ −3/2 µa −5/2 . This property of the resonant islands allows us to automatically detect these regions by computing numerically the derivative of the proper mean motion with respect to the initial semi-major axis a0 . Most of the resonances are then easily identified using the decomposition in continued fraction of n/n 8 . For all the considered cases, this method provided a single rational number − p1 / p2 , thus corresponding to the resonance p1 n 8 + p2 = 0. In Fig. 1, the plateaus that are clearly identified are respectively associated to the resonances (9 : −14), (7 : −11), (5 : −8), (3 : −5), (7 : −12), (4 : −7), (5 : −9), and (1 : −2) with Neptune. Between

FIG. 1. Frequency curve. Frequency ratio n/n 8 versus initial semimajor axis a0 of the test particle. The plateaus correspond to mean motion resonances with Neptune.

8

ROBUTEL AND LASKAR

these plateaus the frequency curve seems to be regular, except for a smaller than 42 AU. These singularities are the signature of chaotic behavior (Section 3). 3. FREQUENCY MAP OF THE SOLAR SYSTEM (i0 = 0◦ )

In order to study the short time dynamics of massless bodies in the Solar System, we consider two different models for the inner system (Mercury to Jupiter) and for the outer planetary system (beyond Jupiter). Particles of inner system are integrated together with all the Solar System’s planets except Pluto over 500,000 years, which corresponds approximately to 2,100,000 orbital revolutions of Mercury and 42,000 for Jupiter. This is largely enough for the FMA to detect the presence of mean motion resonances and to measure the diffusion rate of the proper mean motion. For the outer system, only the four giant planets are taken into account while the masses of the telluric planets are added to the solar one. In this case, the total integration time is 2 Ma, which corresponds to 169,000 orbital revolutions of Jupiter, 12,100 of Neptune, and 2000 of an object at 100 AU. The initial positions and velocities of the planets and Sun are taken from DE200 (Standish 1990) at the date J2000, and expressed in a coordinate system related to the invariable plane of the eight planets. The elliptic elements of the planets and massless bodies are referred to this plane. The initial heliocentric canonical elements (Laskar and Robutel 1995) in this reference frame are given in Table I, together with the measured diffusion rate σ per million years. The initial conditions for the massless particles are taken in the (a, e) plane defined by t0 = J 2000, i 0 = M0 = ω0 = Ä0 = 0. In the (a, e) plane, initial conditions (a0 , e0 ) are chosen on a grid, depending on the location of the particle. The eccentricity e0 is always sampled from 0 to 1, with a constant stepsize of 0.0125, that is, over 80 different values. Then for a0 , 200 values are chosen on a regularly spaced grid in each of the intervals [0.38, 0.9] AU (Fig. 2a), [0.9, 2.0] AU (Fig. 2b), [2, 5] AU (Fig. 2c), and [5, 10] AU (Fig. 2d), and the remaining part is sampled evenly with a stepsize of 0.05 AU for a0 ∈ [10, 100] AU. We have thus numerically integrated altogether 192,000 particles.

The results of all these integrations are plotted in Figs. 2a–2h, where each initial condition corresponds to a colorized dot. For a given particle, three different dynamical situations can occur: (i) The particle is ejected or undergoes a close encounter with a planet before the end of the integration; the corresponding dot is colored in white (the same color as the background). (ii) The particle survives the integration but is in mean motion resonance with a planet (see Section 2.6 for its characterization): a black dot is plotted. (iii) The trajectory is not resonant. The dot is then colored according to its proper mean motion diffusion rate, from deep blue, for motion close to quasiperiodic (σ ≤ 10−8 ), to red for strongly irregular motions (σ ≥ 0). Furthermore, lines of planetary collision are plotted in purple. Their equations are a j (1 + e j ) = a(1 − e) if a ≥ a j and a j (1 − e j ) = a(1 + e) is a ≤ a j , where e, a, e j , a j are respectively the initial eccentricities and the semi-major axis of the particle and the considered planet. Figures 2a–2h provides thus a striking global view of dynamics in the Solar System, which emphasizes the separation of the system into three different regions: The inner Solar System, from 0.38 to 5 AU (Figs. 2a–2c), the outer Solar System, from Jupiter to Neptune (Figs. 2d and 2e), and finally the region beyond Neptune, generally known as the Kuiper belt (Figs. 2f–2h). In the next sections, we will review the two outer regions, for which the present study is the most relevant, and then more briefly the inner planet region. 3.1. Accuracy of the Diffusion Rate and Extent of the Resonances In this study, the FMA is conducted over a relatively short time interval of length T . To be convinced that this time length is sufficient, the stability of the results with respect to T is analyzed by changing the value of T for a selected part of Fig. 2f, namely the rectangle of initial conditions a0 ∈ [45, 50], e0 ∈ [0, 0.5]. For better clarity, the grid in a was refined, with a step size of 0.025 instead of 0.05 in Fig. 2. The results are displayed in Fig. 3. Figure 3a displays the results obtained with the original time

TABLE I Initial Heliocentric Canonical Elements of 8 Planets of the Solar System in the Invariable Plane Reference Frame with Origin J2000 Planet

a (AU)

e

i (◦ )

Ä (◦ )

ω (◦ )

M (◦ )

log σ

Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune

0.38726 0.72387 0.99998 1.52197 5.19109 9.54995 19.21606 30.10826

0.20514 0.00708 0.01668 0.09250 0.04651 0.05360 0.04641 0.00829

6.34230 2.19485 1.57857 1.68041 0.32664 0.93037 1.02832 0.72614

32.48533 51.49619 284.04669 353.16889 313.64397 120.40192 308.62715 191.19602

41.57720 82.85886 175.9095 338.94292 57.34964 328.28190 220.11778 202.46953

174.70400 44.13608 356.98422 19.78240 19.84411 317.81109 140.93673 267.65781

−5.09 −4.84 −5.91 −5.57 −6.53 −6.46 −5.71 −5.58

Note. σ is the diffusion rate measured by million years.

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

9

FIG. 2. Frequency map of the Solar System in the plane: (a0 , e0 ) for i 0 = M0 = ω0 = Ä0 = 0◦ : The color code is for log σ (unit: Ma−1 ). From purple to blue, the motion is regular, from yellow to red it is strongly chaotic.

10

ROBUTEL AND LASKAR

FIG. 3. Frequency maps computed with two different integration times. 2 Ma for Fig. 3a, and 8 Ma for Fig. 3b.

interval with T = 2 Ma, while Fig. 3b was obtained with T = 8 Ma. For a better comparison of these two pictures, the diffusion rate of the particles of the second figure has been rescaled linearly (σ |4 Ma = σ |1 Ma + log(4)), which assumes that the diffusion is linear with respect to time. It is clearly visible that there are barely no changes between the two figures, and that T = 2 Ma was already sufficient for obtaining the major dynamical features of the system. The most visible change is a different measured value for the diffusion rate in the vicinity of the main resonance, which is expected from the asymptotic formula (10). 3.2. Outer Planets Region The outer planets interplanetary region has the characteristic of being vastly depleted: a few particles survive the 2 Ma of integration and almost all the survivors have a very strong proper mean motion diffusion rate, even for small eccentricity (σ > 10−1 ). This instability is essentially due to the overlap of the mean motion resonances caused by the proximity to the giant planets (see Section 3.2.1). In this region, we do not identify many isolated resonances. An isolated resonance island is easily detected by the constancy of the proper mean motion (see Section 2.6), but in a region of strong resonance overlap, the proper mean motion oscillates wildly or jumps between several resonant values. This is why, in Figs. 2d and 2e, only few regions are colored black. The largest of these identified resonant regions corresponds in fact to stable areas surrounding the Lagrangian points L 4 and L 5 of the planets (see also Holman and Wisdom (1993)). The Jupiter’s, Uranus’s, and Neptune’s Lagrangian stable regions are crossed respectively at about 5.2, 19.2, and 30.1 AU. In contrast, close to Saturn’s location, only extremely irregular motions have been detected. This is probably due to the existence of a very unstable region surrounding Saturn’s triangular Lagrangian points (Holman and Wisdom 1993, De la Barre et al. 1996). It is important to point out that the choice of initial conditions for the particles was not designed here for the specific study of the stability of the L 4 and L 5 points. In Fig. 2, the par-

ticles initial mean longitude is equal to 0◦ independently of the longitudes of the planets, but we can deduce from Table I that the relative initial mean longitudes of the bodies (δλ j = λ − λ j ) are respectively equal to δλ5 ≈ −31◦ , δλ6 ≈ −47◦ , δλ7 ≈ 50◦ , and δλ8 ≈ 59◦ . The fact that these quantities are not far from 60◦ (L 4 ) and −60◦ (L 5 ), especially for Neptune, explains why our section of initial conditions intersects the Lagrangian librational regions. The other resonances that are quite well determined are all orbital resonances of order one. The only resonance of order one with Saturn is the (2 : −3), close to 12.5 AU. Three resonances with Uranus have been identified, namely the (6 : −5) at 17.3 AU, the (4 : −5) at 22.3, and the (3 : −4) at 23.2 AU. The last two resonances surrounded by regular areas are resonances with Neptune: the (5 : −4) at 26 AU and the (6 : −5) at 26.6 AU (see Table IV for more details about the location of these resonances). In this large chaotic zone between the giants planets, possibly stable regions are visible. In the first region, between 14 and 14.7 AU with eccentricity lower than 0.1, we found trajectories with diffusion rates lower than 10−4 but always greater than 10−5 . According to Holman and Wisdom (1993) in this region the particles’ lifetime can reach 20 million years while the usual lifetime between Saturn and Uranus is less than 107 years. We find other possible long-time stable candidates at about 25.7 Au and 26.2 AU. The diffusion rate of these bodies is between 10−5.5 and 10−4 . The fact that these two small islands of stability correspond to the “possible long-lived belt of objects between Uranus and Neptune” identified by Holman (1997) shows that the quantity σ , determined by a short time integration (2 Ma), is a good indicator of stability over several hundreds of million years. 3.2.1. Smaller Planetary Masses. In order to illustrate the importance of the mean motion resonances in the interplanetary region, we studied the dynamical behavior of this region for a smaller value of the planetary masses. Figure 4 presents the same results as in Figs. 2d and 2e, but for planetary masses 10 times lower. For this smaller value of the planetary masses, the resonances are narrower. They no longer overlap, and can thus be easily identified. About 30 of these resonances are colored black, and the location of many others is revealed by a narrow vertical line with larger diffusion rate. Each resonance being very close to the next, the enlargement of the width of the resonances by the increase of the perturbation size (this width should increase with the square root of the involved planet’s mass) will lead to its complete overlap, generating the very large scale chaos visible in Figs. 2d and 2e. In contrast to the actual Solar System, for small inclinations, we observe in Fig. 4a a stable region associated with the Lagrangian points of Saturn. 3.3. Kuiper Belt The domain that extends beyond Neptune (the Kuiper belt) is very different from the precedent. Indeed, the dynamics of Kuiper belt objects is essentially governed by the presence of

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

FIG. 4. Frequency map for planetary masses divided by 10. (a0 , e0 ) plane with i 0 = M0 = ω0 = Ä0 = 0◦ . The color code is for log σ (units of Ma−1 ).

Neptune. Figures 2f to 2h show in this region the existence of three different dynamical domains. (i) Above the collision line lies a region that is strongly chaotic, where σ > 10−1.5 (orange to dark red), which we will call the collisional region. (ii) Under a line (in red) parallel to the collision line, starting at e = 0, a = 37 AU, lies a domain of moderately chaotic to regular behaviors where σ < 10−6 is colored green–blue to dark blue, which we will call the regular region. (iii) Between these two regions is an intermediate domain with chaotic trajectories that seems to correspond to the resonance overlap domain, and which we will thus call the resonance overlap region. These three domains with very different dynamical properties are connected and largely penetrated by the mean motion resonance islands (in black), which cover a large area of the intermediate band. Before going further, we will study more precisely these resonances.

11

3.3.1. Kuiper belt resonances. The resonance libration areas have a characteristic V shape in the (a, e) plane: their width is zero for e = 0 while it is maximum in the vicinity of the collision line. When they cross the collision lines, the top part of the resonances has a reversed V shape. A resonance is identified with our procedure when at least two points fell into the libration island. Therefore, as the base of the resonant islands is very narrow, for small eccentricity, we have no points in them, and the black regions in Fig. 2 do not usually reach the axis e = 0. Nevertheless, the presence of the resonances is still detectable by a narrow vertical band where the measured diffusion rate is higher than that of the local background. More precisely, in this case, the measure of the diffusion rate is more a measure of the distortion of the frequency map in the vicinity of a resonance that may be isolated, and thus still a stable region. Many resonances are detected with this method, but without a specific analysis they are often difficult to identify when they are not associated with an identified plateau (in black in Fig. 2). Moreover, the chosen plane of initial condition may not cross any libration zone (elliptic part) of a given orbital resonance (situation discussed in Section 4), but instead crosses the hyperbolic part of the resonance, which will lead also to a local increasing of the measured diffusion rate σ . These phenomena explain the presence of numerous quasivertical lines or narrow bands in the pictures; they also reveal the physical presence of high-order mean motion resonances. At this point, it must be emphasized that the location, the shape, and the width of the resonances are not in this study an analytical or numerical prediction but the effective measures of a section of the resonant areas. These pictures provide a dynamical image of the system in the sense that all initial conditions taken in a black region correspond to resonant trajectories, which may eventually be subject to small diffusion, and thus not necessarily resonant for all time. 3.3.2. Comparisons with analytical computations. It is interesting to compare our determinations of the actual resonance locations in the kuiper belt with those predicted by other means. One of the more complete predictions, obtained by a numerical computation of a resonant normal form of order one, is given in Fig. 1 of Morbidelli et al. (1995). All resonances with Neptune that are present in Morbidelli et al. (1995) are encountered in our study, except perhaps the (5 : −6) close to 34 AU (this resonance is observed in Fig. 9 for i 0 = 30◦ and in Fig. 13a), but this location corresponds to a region of quite nonregular motion that is depleted in less than 107 years (Levison and Duncan 1994). Between these large resonances we detected the presence of many other resonances that are dynamically active, in particular, the (4 : −11) at 59.11 AU between the (3 : −8) and (1 : −3), the (4 : −13) close to 66.15 AU, the (3 : −11) between 71 and 72 AU, and especially on each side of the (2 : −9) located at 82.5 AU the (3 : −13) and (3 : −14) mean motion resonances. Table II gathers all the mean motion resonances well identified with Neptune in this region.

12

ROBUTEL AND LASKAR

TABLE II Mean Motion Resonances with Neptune Clearly Identified in Fig. 2 (i 0 = 0) a0 (AU)

p1

p2

w (AU)

a0 (AU)

p1

p2

w (AU)

30.10 33.40 34.40 34.95 35.55 36.45 37.25 37.70 38.20 38.45 39.50 40.45 40.70 41.20 41.60 42.30 43.15 43.70 44.35 44.55 45.50 45.80 46.05 47.80 49.60 49.80 50.10 50.95 51.30 51.70 52.95 53.65 54.00 54.45 54.70

1 6 9 4 7 3 8 5 7 9 2 9 7 5 8 3 7 4 9 5 7 8 9 1 9 8 7 5 9 4 3 8 5 7 9

−1 −7 −11 −5 −9 −4 −11 −7 −10 −13 −3 −14 −11 −8 −13 −5 −12 −7 −16 −9 −13 −15 −17 −2 −19 −17 −15 −11 −20 −9 −7 −19 −12 −17 −22

0.70 0.10 0.10 0.50 0.25 0.75 0.20 0.20 0.25 0.15 0.95 0.10 0.15 0.25 0.20 0.80 0.10 0.50 0.25 0.30 0.05 0.15 0.10 1.65 0.15 0.15 0.15 0.25 0.15 0.55 0.75 0.20 0.25 0.10 0.20

55.45 56.60 57.00 57.35 57.90 59.15 59.60 59.90 60.70 60.95 61.15 62.60 64.50 65.50 66.15 67.25 67.90 68.20 68.80 69.45 70.85 71.65 72.80 73.45 75.90 78.55 79.10 80.10 81.05 82.15 84.20 85.25 85.90 88.10

2 7 5 8 3 4 9 5 7 8 9 1 8 5 4 3 8 5 9 2 5 3 4 5 1 5 4 3 5 2 3 4 5 1

−5 −18 −13 −21 −8 −11 −25 −14 −20 −23 −26 −3 −25 −16 −13 −10 −27 −17 −31 −7 −18 −11 −15 −19 −4 −21 −17 −13 −22 −9 −14 −19 −24 −5

1.00 0.15 0.20 0.15 0.80 0.50 0.15 0.15 0.20 0.10 0.10 1.60 0.10 0.20 0.45 0.75 0.15 0.15 0.10 0.95 0.15 0.75 0.45 0.25 1.80 0.25 0.45 0.80 0.15 1.05 0.80 0.45 0.20 1.65

Note. a0 corresponds to the location of the base of the libration island. The integers p j define the resonance, and w is its approximate maximal width measured along a horizontal line. Here the measured width w of a resonance depends on the chosen section. Thus, w gives only a lower bound of the maximal width of the resonances.

If we compare the shapes of the libration areas in these two different models,we find that they are very similar except for the symmetry with respect to a vertical line and the width. Indeed the libration areas in Morbidelli et al. (1995) are almost symmetric, which is not the case in our figures. This asymmetery increases with distance from Neptune, which is particularly well illustrated by the (1 : −5) resonance for which most of the libration island is on the right side of the vertical line a = 88.1 AU (see also Fig. 12f). This phenomenon will be discussed in more detail in Section 4. On the other hand, the maximal width (measured along a horizontal line) of the libration areas (Tables II and IV) is usually smaller in our study than in Morbidelli et al. (1995). Indeed, Morbidelli et al. (1995) provided an upper bound of the domain of libration, considering its projection in the plane (a, e)

independently of the angles and not a section of this domain. It is remarkable that even far away from Neptune, the mean motion resonances occupy an important region of the phases space and play a nonnegligible dynamical role. This could allow for example trapping of particles in resonances at large eccentricity under some dissipative effect (collisions or planet migration). In this region, even for large values of the semi-major axis, mean motion resonances with Uranus are observable. Indeed the (1 : p2 ) resonances with Uranus are recognizable by a distortion of the frequency map (small vertical band of high diffusion) close to the right edge of the (2 : p2 ) resonance with Neptune. The (1 : −2) is marked by a few black points near the tringular Lagrangian stable region of Neptune, the (1 : −3) is located on the right side of the (2 : −3) resonance with Neptune (see Section 6 and Table II), and so on up to the (1 : −10) resonance at 89.2 AU. In the next section, we will study in more detail the three different dynamical domains of the Kuiper belt. 3.4. Kupier Belt Dynamical Regions 3.4.1. The collisional region. The main characterstic of the domain above the collision line with Neptune is its strongly irregular dynamics. The particles that have not being ejected over a 2-Ma timescale have a very high proper mean motion diffusion rate, greater than 10−1.5 for almost all particles. This strong chaos is mostly due to close encounters with Neptune, and secondarily with the other planets. Almost all particles whose initial elliptic orbits cross Saturn’s orbit are ejected before 2 Ma. On a longer time span, these bodies will probably be removed from the Solar System after planetary close approach. The only exception to this situation is, as for Pluto, when the particle evolves in the stable area associated with a orbital resonance with Neptune. This situation provides a mechanism that prevents the particle from collision with the considered planet, at least when the amplitude of libration of the critical angle is sufficiently small. This is why many resonances penetrate the large chaotic zone above the collision line. As many as 52,751 particles were initially in the collision domain. We found that 40.5% of these particles were removed from the integration before 2 Ma. The number N (t) of removed particles versus time is plotted in Fig. 5. This very regular curve is perfectly well approximated by the function N (t) = 9888 log t − 41014 for t which was fitted to the data over the interval [4 × 105 , 2 × 106 ] years. Thus, as in Holman and Wisdom (1993), we find that the number of Neptune encounters decays as 1/t (see Sim´o et al. (1995) for a rigorous estimate of the escape time). Before this instant, the decay is faster: the particles that are initially close to the giant planets are rapidly eliminated. The clearing of this region probably continues at the same rate until all particles that are not in resonance with Neptune have been removed. Particles in stable libration resonant island should survive, and particularly at the tringular Lagrangian points with Neptune. Nevertheless, these regions are probably also slowly erobed by secular dynamics (Morbidelli et al. 1995, Grazier et al. 1999b).

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

13

between 40 and 42 AU) is the domain where secular resonances interact (Morbidelli et al. 1995). This leads to a possible long timescale chaos that would reduce the particles lifetime (Holman and Wisdom 1993, Levison and Duncan 1993). For a semi-major axis greater than 45 AU, the secular frequencies of the massless body are too small (at least outside mean motion resonances) to lead to low-order secular resonances with the giant planets. For this reason, we can suppose that the short-time stable regions, having a diffusion rate lower than 10−6 or 10−7 , are stable for a very long time.

3.5. Inner Planets Region

FIG. 5. Number of ejected bodies versus integration time in million years (top), and difference with the fitted model N (t) = 9888 log t − 41014 (bottom).

3.4.2. Resonance overlap region. The intermediate region is perhaps most interesting. This domain, colored yellow to green–yellow, corresponds to a mean motion diffusion rate between 10−2.5 and 10−1.5 . Figure 2 shows clearly that mean motion resonances connect in this band. Here the diffusion, which is smaller than in the above one, is essentially driven by resonance overlap. The chaotic motions induced by this phenomenon will probably lead to a complete depletion of this area by close encounters. For the collision domain, this region will certainly be cleared for sufficiently long times. If we look carefully to the overlap domain, we can see that it is bounded by two parallel curves: the upper bound is given by the collision line with Neptune, while the lower one is a curve of equation a(1 − e − δ) = q0 , where q0 is the initial perihelic distance of Neptune and δ a parameter that can be estimated in Fig. 2 as δ ≈ 0.196. This observation is in agreement with the overlap criterion of Wisdom (1980), which was developed for the planar circular restricted three-body problem. In the case of Neptune, the application of Wisdom’s criterion would predict that the region of resonance overlap lies from 27.4 to 32.8 AU, which slightly underestimated its size. Indeed the observed chaotic region (Fig. 2) ranges from about 27 AU to about 33 to 35 AU. This small underestimation can be explained as our model contains four planets whose orbits are neither circular nor planar. The agreement is still acceptable considering also that the resonances taken into account by Wisdom for the establishment of his criterion were only of order one. 3.4.3. Regular region. The last dynamical domain pointed out in our study is the one of slow diffusion. This domain in blue contains the orbits with σ < 10−6 . For a < 45 AU, this region is limited to very small spots with low eccentricity. These areas are not connected, and the stability of these orbits may need further studies. Indeed, this region (close to 36 AU and

As this region has already been studied in detail, in particular for the dynamics of asteroids in mean motion resonances (see Nesvorny and Ferraz-Mello (1997) and references therein), we will limit ourselves to a brief discussion of the results in the case of the inner Solar System (Figs. 2a to 2c). In this study, as the time of integration is only T = 0.5 Ma, the diffusion rate σ is rescaled by addition of the factor log 4 (Section 3.1), and as for the outer system, it is given per million years. As for the outer system, the collision regions are extremely unstable. Between Mars and Jupiter, we recognize the resonant structure of the asteroid belt, and all the principal mean motion resonances with Jupiter are identified by their libration region (in black) for (2 : −1), (5 : −3), and (3 : −2), or by a vertical band filled with very unstable or ejected points, especially for (9 : −2), (4 : −1), (3 : −1), (5 : −2), (7 : −3), and (4 : −3). In this region, the stepsize in a being equal to 0.015 AU is too small to detect the libration islands of very narrow resonances, particularly the three-body mean motion resonances whose typical width is at most a few 10−3 AU (Nesvorny and Morbidelli 1998). Even if these thin structures are not directly identified by our automated procedure, some are present in the pictures as narrow vertical lines. In this inner region, only a few number of resonances with terrestrial planets have been detected. If we assume that the width of a resonances is roughly proportional to the quanp tity m j /m 0 where m 0 is the Solar mass, the expected width of the resonances with Venus or the Earth is approximately 20 times smaller than with Jupiter, and 60 times for Mercury and Mars. If we compare with the (2 : −1) resonance with Jupiter, we can estimate that the width of the same resonance with terrestrial planets is smaller than 0.005 AU, a quantity which is on the order of the resolution in the semi-major axes for Figs. 2b and 2b. Below the collision lines, only a very few areas seem to be regular. Except for well-known regions of the asteroid belt, only a small region close to 0.6 AU and that between 1.2 and 1.4 AU could be stable for a long time. Evans and Tabachnik (1999) have shown that in this region, bodies that are on circular orbits may survive for a long time. This region certainly needs to be studied further, with higher resolution, in order to obtain a clearer view its dynamics. The benefit gained from the present study is already for this region to rule out a large part of the phase space where no regular orbits can exist.

14

ROBUTEL AND LASKAR

4. GEOMETRY OF THE MEAN MOTION RESONANCES IN THE PHASE SPACE

In the previous section, the resonant areas have been studied in a (a0 , e0 ) plane. Now, we will study the effect of initial phases variation on the frequency map. In order to understand this effect, we have concentrated our attention to the phase space region where the semi-major axis is contained between 40 and 50 AU. We fixed the eccentricity value to e0 = 0.25 and let one of the phase angles vary (ω0 ∈ [0, 360] and M0 = 0 in Fig. 6a, M0 ∈ [0, 360] and ω0 = 0 in Fig. 6b.), and integrated 18,000 particles over 2 Ma. The stepsize is 0.05 AU for a0 and 4.5◦ for the angles. The bottom lines of the two figures (ω0 = 0 in 6a and M0 = 0 in 6b) are the same as the the segment e0 = 0.25, a0 ∈ [40, 50] AU of Fig. 2f. As this line is crossing the collision curve at 41 AU, it explains why the regularity of the dynamics increases from left to right in Figs. 6a and 6b. These figures, which represent sections of the phase space, clearly show the phase dependence of the libration islands associated with resonances. We recognize very easily in Figs. 6a and 6b vertical chains of elliptical islands separated by regions of hyperbolic dynamics. This chains are distorted by the vicinity of low-order resonances. This phenomenon is particularly rele-

vant in the neighborhood of the (1 : −2) resonance (a0 close to 48 AU). 4.1. Simple Model of Resonance These figures can be interpreted using a simple model of resonance. As the resonances identified in Fig. 6 are all mean motion resonances with Neptune, the simplest model that we can consider is the circular restricted three-body problem. Its Hamiltonian can be written as µ2 + n 8 38 232 X + A(k1 ,k2 ) (3, 38 , x, x¯ , y, y¯ ) exp i(k1 λ8 + k2 λ) (16)

H=−

k∈Z 2

√ √ p √ √ 2 exp(i$ ), y = with 3 = µa, x = 3 1 − 1 − e 3 q √ 1 − 1 − e2 (1 − cos i) exp(iÄ), λ8 = n 8 t + λ0 , and 38 is the conjugated variable of 38 . The coefficients A(k1 ,k2 ) are sum of monomials of the form α α¯ β β ¯ y y¯ , Ck1 ,k2 ,α,α,β, ¯ β¯ x x ¯

(17)

where the integers α, α, ¯ β, β¯ satisfy the Dalembert relation k1 + k2 = α¯ − α + β¯ − β (see Laskar and Robutel (1995) for details). As r √ 3 x≈ (18) e, y ≈ 23 sin(i/2) 2 the coefficient (17) is equivalent to K

m 8 (α+α) ¯ e ¯ sin(i/2)(β+β) , m0

(19)

where K is a constant. The Dalembert relations imply that for quasi-planar motions (which is the case in this section), the resonant angles associated with a ( p1 : p2 ) mean motion resonance with the jth planet (angles that remain in the resonant normal form) correspond to the harmonics of the combination p1 λ8 + p2 λ − ( p1 + p2 )$.

(20)

In contrast, for bodies orbiting far from the invariable plane and in quasi-circular motion (this situation will be encountered in Section 6), the resonant angles are harmonics of p1 λ8 + p2 λ − ( p1 + p2 )Ä

FIG. 6. Effect of phases variation on the geometry of the mean motion resonances. Initial conditions of the test particles are e0 = 0.25, i 0 = 0◦ , Ä0 = 0◦ for the three figures, with (a0 , ω0 , M0 ) ∈ [40, 50] × [0, 360] × {0} for Fig. 6a. (a0 , M0 , ω0 ) ∈ [40, 50] × [0, 360] × {0} for Fig. 6b, (M0 , ω0 , a0 ) ∈ [0, 360] × [0, 360] × {42.2} for Fig. 6c.

(21)

and in the general case, when neither eccentricity nor inclination is small, the resonant angles correspond to the harmonics of the two-parameter family kp1 λ8 + kp2 λ + (α − α)$ ¯ + (α¯ − α − kp1 − kp2 )Ä, (22) where k and α − α¯ are the two integer parameters.

15

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

In this section, as the particles are initially located in the invariable plane, the planar-restricted problem provides a sufficient approximation to give a qualitative interpretation of Figs. 6a and 6b. If we focus on the ( p1 : p2 ) mean motion resonance with Neptune, a classical way to proceed is to perform the time-dependent canonical transformation: ϕ1 = p1 λ8 + p2 λ − ( p1 + p2 )$ ϕ2 = α1 λ8 + α2 λ − (α1 + α2 )$ F1 = (α1 + α2 )3 − a2 x x¯

(23)

F2 = −( p1 + p2 )3 − p2 x x¯ with FIG. 7. Phase portrait of the Hamiltonian H = J 2 /2 + cos φ + cos 2φ. Compare the shape of the island to the island of the 1/2 resonance in Figs. 6a and 6b.

p1 α2 − p2 α1 = 1, which after an averaging of first order on the angle ϕ2 leads to an approximation of the normal form of the previous Hamiltonian µ + n 8 ( p1 F1 + α1 F2 ) 2( p2 F1 + α2 F2 )2 X A0k (F1 , F2 ) cos(kϕ1 ). +

would correspond to the hyperbolic fixed point of (26) associated with the ( p1 : p2 ) resonance, and (A01 is a negative constant)

2

H0 = −

(24)

k≥0

4.2. Location of the Islands The number of harmonics taken into account in the previous sum depends on the chosen resonant model (from one in the pendulum approximation to an infinity for more accurate models) (Poincar´e 1902, Henrard and Lemaitre 1983, Robutel 1993, Morbidelli et al. 1995). For a qualitative understanding of the figures, we will use the simple pendulum model, obtained by retaining only the first harmonic (it is assumed to be the term of largest amplitude, which is true for sufficiently small eccentricity) of (24), and expanding H0 in the neighborhood of the resonance up to order 2 for the first two terms of (24) and zero for the third. If we put F1 = F1∗ + I , and F2 = F2∗ where F2∗ is the value of the first integral F2 and F1∗ the value of F1 at the resonant, that is, µ2 p 2 + n 8 p1 = 0, with 3∗ = p2 F1∗ + α2 F2∗ , 3∗3

(25)

3 p22 µ2 2 I + A01 (F1∗ , F2∗ ) cos(ϕ1 ). 2 3∗2

(26)

If we translate the well-known dynamics of this model in the initial variables, we found that the initial values ¡

(a, ϕ1 ) = ( p2 / p1 )

2/3

¢

a8 , 2kπ , k ∈ Z

(27)

(28)

to the elliptical one. With more sophisticated models, others fixed points may appear. In particular, for large eccentricities, bifurcations can arise, changing the nature of the elliptical point at ϕ1 = π and leading to the birth of two new elliptical fixed points in the libration island (Message 1958, Beaug´e 1994). This is clearly visible is the phase portrait of the simple Hamiltonian (Fig. 7) H=

J2 + cos φ + cos 2φ. 2

(29)

In fact, in the complete, five-body restricted problem, the points that can be found in Figs. 6a and 6b do not correspond to fixed points but to initial conditions of quasiperiodic orbit lying on a resonant torus of lower dimension than the usual KAM tori of quasiperiodic orbits. Nevertheless the pendulum approximation in the restricted circular problem is sufficient to understand the location of the islands in Fig. 6. From (20) we can express the critical angle by ϕ1 = p1 λ8 + p2 M − p1 (ω + Ä).

we get, after expansion the pendulum Hamiltonian, Hpend = −

¡ ¢ (a, ϕ1 ) = ( p2 / p1 )2/3 a8 , (2k + 1)π , k ∈ Z

(30)

As the initial longitude of the node Ä and the mean anomalies M of trajectories plotted in Fig. 6a are equal to 0, from (30) the fixed points of the pendulum corresponding to the ( p1 : p2 ) resonance are located close to the p1 points a = ( p2 / p1 )2/3 a8 , ω = λ8 −

2kπ p1

(31)

16

ROBUTEL AND LASKAR

TABLE III Computed Prediction for the Hyperbolic Points of Fig. 5 for the Selected Resonances (1 : −2), (3 : −5), and (4 : −7) p1

p2

a

ωh

ωe

Mh

Me

1 —

−2 —

47.791 —

301

121

150 330

60 240

3 — — — —

−5 — — — —

42.32 — — — —

61 181 301

1 121 241

37 109 181 253 325

1 73 145 217 289

4 — — — — — —

−7 — — — — — —

43.7 — — — — — —

31 121 211 301

76 166 256 346

18 69 121 172 223 275 326

43 95 146 198 249 301 352

Note. The coordinates (a, ωh,e ) correspond to Fig. 6a and (a, Mh,e ) to Fig. 6b (the index h corresponds to the hyperbolic points and e to the elliptic ones).

for the hyperbolic ones, and close to the p1 points a = ( p2 / p1 )2/3 a8 , ω = λ8 −

(2k + 1)π p1

(32)

for the elliptic one. Identically, as in Fig. 6b where the quantities ω and Ä are fixed to 0, coordinates of the 2 p2 fixed points of the model pendulum are respectively 2kπ p1 − λ8 p2 p2

(33)

2(k + 1)π p1 − λ8 p2 p2

(34)

a = ( p2 / p1 )2/3 a8 , M = for the hyperbolic ones, and a = ( p2 / p1 )2/3 a8 , M =

In Fig. 6c, we play another game. The initial action-like variables are fixed (a0 = 42.2 AU, e0 = 0.25, and i 0 = 0), the initial longitude of the node is fixed (Ä0 = 0), while the two others angles M and ω are varied in the interval [0, 360] with a step size of 4.5◦ . As expected, the resonant area, corresponding to the (3 : −5) orbital resonance with Neptune, is confined in a band of slope −5/3. Indeed, in the libration region, the critical angle ϕ1 satisfies the inequalities |ϕ1 − π| ≤ A(a, e),

(35)

where the libration amplitude A, which depends on a and e, is approximately equal to 2π/3. Then the libration band upper and lower bounds are given by the relations 5M + 3ω = π + 3λ8 ± A. From expression (31), we can conclude that the line M0 = Ä0 = i 0 = 0, e0 = 0.25, ω0 = λ8 parametrized by a0 , crosses only hyperbolic regions. This phenomenon can be observed in Fig. 8b where the line ω0 = λ8 ≈ 301◦ does not cross resonant islands. In order to illustrate this phenomenon, we present in Fig. 8b the same map as in Fig. 2, but with initial argument of perihelion ω0 = 301◦ instead of ω0 = 0. This figure is comparable to Fig. 8b, which is a magnification of the rectangle (a0 , e0 ) ∈ [40, 50] × [0, 0.5] of Fig. 2f. The global aspect of these two maps is essentially the same: we find the same extremely chaotic behavior above the collision line with Neptune. Just under the collision line lies the resonance overlap region, and beyond about 45 AU, the regular region. In fact, the only differences are in the width and in the dynamical nature of the crossed resonant areas: The three libration regions associated with the (3 : −5), (4 : −7), and (1 : −2) mean motion resonance with Neptune almost vanish for ω0 = 301◦ . In counterpart, their associated unstable region is always present. Moreover, the areas filled with trajectories of small diffusion (σ = 10−5 ) have approximately the same shape and the same size in the two figures, even close to resonances of small order. It is in this sense

for the elliptic ones. To illustrate the previous discussion, the predicted location of the points corresponding to the pendulum fixed points (formulas (31) to (34)) are given in Table III for the resonances (1 : −2), (3 : −5), and (4 : −7). 4.3. Size of the Islands In counterpart, the pendulum approximation is not able to predict accurately the width of the libration areas. Indeed, the shape of the (1 : −2) mean motion resonance with Neptune is not that expected from a simple pendulum model. Figures 6a and 6b show clearly that the amplitude of the second resonant harmonic (2ϕ1 ) is on the same order as the first one. One can indeed see that the shape of the island in Figs.6a and 6b is very similar to that in Fig. 7, which thus seems to be sufficient to reveal the bifurcation of the center of the island into an hyperbolic orbit (see Morbidelli et al. (1995)).

FIG. 8. Effect of the phase variation on the width of the resonances in the (a0 , e0 ) plane. (a) ω0 = 0◦ . (b) ω0 = 301◦ .

17

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

that we can say that the results obtained using FMA, except inside the libration regions, are independent of the initial phases. It is precisely this observation (see Laskar (1999)), which makes it possible to obtain a complete view of the dynamics of the phase space, despite the fact that we have fixed all the phase angles.

5. FREQUENCY MAP FOR HIGH INCLINATION

In Section 3 the test particles were initially in the invariable plane of the giant planets (i 0 = 0). Now we will perform the same study, but for an initial inclination i 0 of 30◦ . The global results are plotted in Fig. 9, which differs only from Fig. 2 by the choice of initial inclination of the particles i 0 . Comparison between these two figures shows several differences: (i) Particles survive a longer time in the collisional region at high inclination. This is understandable as for highly inclined orbits, close encounters with planets are less frequent than with low inclination. Thus, we must wait a longer time to see the effects of instabilities generated by close approaches, but after sufficient time, all the bodies above the collision lines that are not inside a resonance libration island will probably be ejected. Indeed Gladman and Duncan (1990), integrating two samples of particles with inclination of 0◦ and 10◦ respectively, have noticed that the depletion time does not really depend on the inclination. (ii) Particles are more chaotic when their initial orbit is just on the collision line. This is particularly visible for the collision line with Neptune. This is probably because these orbits will be much more perturbed than the orbits than are above the curve, which will have higher relative velocity during planetary encounters and will thus be less affected. (iii) Diffusion rate is in general higher below the collision lines than in the planar problem. One reason is the increase of amplitude of the size of many mean motion resonances that were not significant in the planar problem. As was stated in Section 4, we have now a complete two-parameter family of resonances. The overlapping of resonances will thus arise for a lower value of the eccentricity of the particles. A second reason, at least for a semi-major axis smaller than 60 AU, is the probable development of secular resonances. Indeed, in case of chaotic behavior generated by secular resonances, we still expect to observe some diffusion of the mean motion frequencies (see Laskar (1993)). We have not attempted here to identify these possible resonances as the time span we used in the present work is not sufficient to obtain precise results on secular resonances, but we expect to do so in further work. Nevertheless, it is already understandable that increasing the inclination will increase the angular momentum deficit (AMD) of the particle (Laskar 1997, 2000), and thus its perturbations in the secular system. (iv) Another interesting phenomenon is the presence of many libration islands associated with mean motion resonances with Uranus, which are clearly seen in Figs. 9f–9h by black regions

TABLE IV Mean Motion Resonances Clearly Identified in Fig. 9 (i 0 = 30◦ ) a0 (AU)

Planet

p1

p2

w (AU)

11.60 12.45 13.40 16.55 17.30 19.15 21.65 22.30 23.20 25.45 25.90 26.65 27.55 29.90 32.55 33.05 33.85 34.00 34.55 35.30 37.10 37.80 40.10 43.00 44.40 45.90 46.50 50.55 51.35 52.65 53.95 56.45 59.05 60.15 61.45 61.90 65.65 69.05 70.90 83.90

6 6 6 7 7 7 7 7 7 8 8 8 8 8 7 7 7 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

3 2 3 5 7 1 5 4 3 9 5 6 8 1 5 4 3 5 5 2 3 4 1 3 2 3 4 4 3 2 3 1 3 2 3 4 4 4 1 1

−4 −3 −5 −4 −6 −1 −6 −5 −4 −7 −4 −5 −7 −1 −11 −9 −7 −6 −12 −5 −8 −11 −3 −10 −7 −11 −15 −17 −13 −9 −14 −5 −16 −11 −17 −23 −25 −27 −7 −9

0.10 0.25 0.25 0.10 0.10 0.55 0.25 0.20 0.20 0.15 0.20 0.20 0.20 0.10 0.15 0.16 0.17 0.15 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39

Note. The resonances with Neptune given in Table II are not reported here. In the second column is the index of the involved planet. a0 corresponds to the location of the base of the libration island. The integers p j define the resonance, and w is its approximate maximal width measured along a horizontal line.

on the line of collision with Neptune (Table IV). In fact, as these bodies are not protected from close encounters with Neptune they will are only be temporarily trapped in these resonances. (v) Let us note finally note that the triangular Lagrangian regions of Saturn is visible, in Fig. 9d, as two vertical bands corresponding to the stable region (the neighborhood of these points is still unstable), but, as mentioned in G´omez et al. (1993), and De la Barre et al. (1996), the stable area surrounding the equilibrium Lagrangian point seems to be larger to high inclination.

18

ROBUTEL AND LASKAR

FIG. 9. Frequency map of the Solar System in the plane (a0 , e0 ) for i 0 = 30◦ and M0 = ω0 = Ä0 = 0◦ . The color code is for log σ (units of Ma−1 ).

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

6. FREQUENCY MAPS OF KUIPER BELT OBJECTS

In this section we study the dynamics of 63 known Kuiper Belt’s objects including Pluto. The selected objects are those from the Marsden’s Kuiper Belt Objects table (Web page: http:// cfa-www.harvard.edu/cfa/TNOs.html), observed at least at two different oppositions. For each object, we make a frequency map in the plane (a0 , ω0 ). The results are plotted in Figs. (10–13), where each small figure provides a real dynamical ID of the considered particle. Table V lists the parameters of the objects: name (first column), heliocentric canonical elliptical elements at the date J2000 (columns 2 to 7), and measured diffusion rate σ per million years (column 8). When the object is in mean motion resonance with Neptune, the integer coefficients of the resonance ( p1 , p2 ) are given in column 9, and in the last column is the index of the corresponding figure. In each figure (10 to 13) the initial elements of the test particles are fixed to those of the considered Kuiper belt object, except for the argument of the perihelion ω0 , which varies from 0◦ to 360◦ with a step size of 4.5◦ , and the semi-major axis a0 for which the range of variation depends of the location of the object. This study will allow us the explore the dynamical vicinity of all the observed objects. It is important to note that knowledge of the actual location and of the dynamical influence of the neighboring resonances will allow us to have a precise idea of the dynamical behavior of the considered body, even if its elliptical elements are still not accurately known. 6.1. Objects in the (2 : −3) Resonance In Fig. 10 are gathered the objects that are in the (2 : −3) resonance with Neptune. In all these figures, a ranges in the interval [38.5, 40.5] with a step size of 0.025 AU. These bodies, except Pluto (Fig. 10u), are classified by increasing initial eccentricity from 0.1 (Fig. 10a) to 0.33 (Fig. 10t). In all cases, the libration zone associated to the (2 : −3) mean motion resonance with Neptune covers a significant area. This resonant region is always composed of two lobes, which shows that the amplitude of the critical angle 2λ8 − 3λ + $ = 2(λ8 − Ä) − 3M − 2ω is dominant (see expressions (20), (21), (22)), even when the bodies are far above the invariable plane (the initial inclination of 1995QZ9 is greater than 19◦ ). The 21 objects of Fig. 10, represented by a white disk, are all in the libration area of the (2 : −3) resonance. For this reason, the diffusion rate σ is not reported in Table V, (see Section 2.6). Some of the Kuiper objects, like 1996TQ66 (Fig. 10c), 1998US43 (Fig. 10d), 1998HK151 (Fig. 101), 1995HM5 (Fig. 10n), or 1996TP66 (Fig. 10t), are deep inside the libration area and therefore will probably remain there a long time, whereas others, like 1998VG44 (Fig. 10m), 1998HH151 (Fig. 10f), and 1995QY9 (Fig. 10p), are close to the edge exterior of this zone where. If the orbital elements of these particles are well determined, this would mean that they are only transitory trapped, and could get lost after a long time, but it could also result from uncertainty in the determination of the orbital elements.

19

From Figs. 10a to 10e the initial eccentricity of the Kuiper belt body is sufficiently small and the particle remains below the collision line with Neptune. Thus, dynamical structures can be observed on a mildly chaotic background. In particular, on the right side of the (2 : −3) resonances chains of resonant islands are distorted by the presence of the major resonance. In the first three figures, one can recognize the (1 : −3) resonance with Uranus, which seems to be composed of at least two libration lobes, but in fact only one lobe is present, which is very much distorted by the presence of the main (2 : −3) resonance. With the initial location of the other bodies represented in Fig. 10 being close or above the Neptune’s collision curve (see Fig. 2), only the (2 : −3) resonance is visible, while all the other details are erased by the great chaotic area surrounding this resonance. 6.2. Kuiper Objects with 41 ≤ a ≤ 45 AU A neighborhood of the 29 bodies with semi-major axis between 41.3 and 44.6 AU is presented in Figs. 11a–11u and 12a– 12h. The semi-major axis a0 spans the interval [41, 45] with a step size of 0.025 AU, while the angle ω0 is the interval [0, 360] with a step size of 4.5◦ . The initial eccentricity of the Kuiper belt’s objects ranges from 0.006 for 1996TK66 (Fig. 11a) to 0.218 for 1994JS (Fig. 12h). The main resonances encountered in this region are (5 : −8) at 41.2 AU, (8 : −13) at 41.6 AU, (3 : −5) at 42.3 AU, (7 : −12) at 43.5 AU, (4 : −7) at 43.7 AU, and (5 : −9) near 44.6 AU. The width of the libration islands of these resonances is clearly related to the size of the initial eccentricity. Except in a few cases, the size of the islands increases from Figs.11a to 12h. Indeed, this phenomenon is also modulated by the initial inclination of the bodies which explains these variations of resonance amplitude: the eccentricity of 1997RX9 and 1998HM151 (Figs. 11h and 11i) are almost the same while their inclination is 29.21◦ for the first body and 1.27◦ for the second. In this region, the main resonance is the (3 : −5) resonance with Neptune (at 42.3 AU), which presents a strange behavior: sometimes three lobes are present, sometimes five, and in other cases neither three nor five. This behavior can indeed be understood when related to the analysis of Section 4. When the eccentricity is dominant compared to the inclination, as is the case for 1998HP151 (Fig. 11u), 1998KR65 (Fig. 12a), 1996TS66 (Fig. 12g), and especially 1994JS (Fig. 12h), the resonance composed of three islands as the amplitude of the critical angle 3λ8 − 5λ + 2$ = 3(λ8 − Ä) − 5M − 3ω is dominant. On the other hand, if the inclination is dominant, five islands are found as for 1998FS144 (Fig. 11c), 1997QH4 (Fig. 11f), and 1997RX9 (Fig. 11h). In this case the amplitude of the angle given in (21) is dominant, that is, the combination 3λ8 − 5λ + 2Ä = 3(λ8 − Ä) − 5(M + ω). The intermediate situation arises when none of the critical angles (see formula 22) is dominant, which is the case for 1998HO151 (Fig. 12d) and 1996TO66 (Fig. 12f) for which the three lobes are very deformed, or for 1995SM55 (Fig. 12c) where four lobes already exist, and for 1996RQ20 (Fig. 12e) where a fifth lobe is almost

20

ROBUTEL AND LASKAR

TABLE V Pluto and Kuiper Belt Objects Retained in Section 6: Initial Conditions, Diffusion Rate, Resonance, and Index of the Corresponding Figure Name

a (AU)

e

i

M (◦ )

ω (◦ )

Ä (◦ )

Pluto 1992 QB1 1993 FW 1993 RO 1993 SB 1993 SC 1994 ES2 1994 EV3 1994 GV9 1994 JQ1 1994 JR1 1994 JS 1994 TB 1994 VK8 1995 DA2 1995 DB2 1995 DC2 1995 HM5 1995 QY9 1995 QZ9 1995 SM55 1995 WY2 1995 YY3 1996 KV1 1996 RQ20 1996 RR20 1996 SZ4 1996 TK66 1996 TL66 1996 TO66 1996 TP66 1996 TQ66 1996 TS66 1997 CQ29 1997 CR29 1997 CS29 1997 CT29 1997 CU29 1997 CV29 1997 QH4 1997 QJ4 1997 RT5 1997 RX9 1997 SZ10 1998 FS144 1998 HH151 1998 HJ151 1998 HK151 1998 HL151 1998 HM151 1998 HN151 1998 HO151 1998 HP151 1998 HQ151 1998 KR65 1998 KS65 1998 SM165 1998 US43 1998 VG44 1998 WA25 1998 WH24 1998 WV24 1998 WX24

39.5552 44.0393 43.7510 39.3477 39.4324 39.6372 45.8118 42.9205 43.7554 44.1243 39.4526 42.3056 39.5992 42.7945 36.3787 46.4089 44.0907 39.5700 39.8223 39.5544 41.7994 46.4899 39.3530 45.2998 43.9570 39.7388 39.6834 42.6804 83.8959 43.3817 39.5256 39.4978 43.9714 45.1260 47.0083 43.8660 43.7220 43.4404 42.1727 42.7514 39.3932 41.7564 41.4172 48.2147 41.8813 39.0207 43.3399 39.4623 40.7833 44.5776 37.9647 41.4504 44.0043 39.4582 44.0832 43.8880 49.3264 39.4354 39.1614 42.5791 45.9145 39.3728 43.3311

.2503 0.072 0.049 0.200 0.320 0.188 0.115 0.041 0.061 0.049 0.119 0.218 0.318 0.030 0.074 0.134 0.063 0.253 0.266 0.148 0.097 0.129 0.219 0.111 0.102 0.183 0.259 0.006 0.583 0.118 0.333 0.124 0.128 0.116 0.213 0.007 0.032 0.039 0.046 0.031 0.220 0.082 0.034 0.369 0.010 0.167 0.040 0.230 0.091 0.037 0.052 0.098 0.078 0.285 0.081 0.059 0.407 0.142 0.246 0.013 0.121 0.108 0.042

15.57 3.09 7.67 3.30 2.91 5.91 1.16 2.28 1.48 3.90 2.71 13.13 13.50 0.94 5.14 2.68 1.66 4.75 5.87 19.35 27.03 0.83 1.19 6.58 31.89 4.62 5.01 2.95 24.55 28.01 7.12 14.90 8.88 1.61 17.72 3.74 0.92 2.63 6.48 13.80 17.39 11.92 29.21 12.15 12.78 8.84 2.03 5.30 28.10 1.27 24.44 15.19 1.34 12.86 0.40 0.40 13.22 11.42 1.61 0.86 11.27 1.90 1.17

14.76 7.06 328.36 3.51 326.85 41.58 283.09 173.10 53.07 305.96 5.30 329.55 332.62 238.49 30.29 3.75 251.17 329.21 351.07 35.37 309.26 269.60 7.10 339.33 5.11 114.31 344.32 225.97 359.35 114.51 359.57 2.49 340.81 39.92 49.82 330.04 215.67 221.52 350.08 358.01 314.92 97.38 240.52 12.23 87.36 347.37 37.97 2.92 43.71 50.13 164.79 211.72 296.03 9.85 91.35 312.49 26.60 37.59 340.61 4.34 315.06 292.44 185.82

113.56 31.72 31.64 163.06 108.79 331.65 11.26 55.92 216.58 273.53 82.00 242.42 102.06 196.92 328.51 350.42 76.46 41.08 37.44 138.23 71.68 176.52 130.97 195.33 1.31 34.96 47.95 127.07 180.98 244.52 81.36 24.40 139.26 312.77 298.87 210.70 319.37 299.29 27.72 11.20 85.00 70.07 303.01 347.98 204.74 12.50 181.72 193.31 134.92 243.58 342.01 214.17 315.72 341.08 259.46 55.09 121.30 130.96 302.44 153.45 70.88 247.72 267.77

107.03 326.82 195.97 191.95 321.23 336.89 240.16 332.15 261.92 358.50 161.69 47.55 310.60 351.98 129.92 137.19 194.14 202.15 325.95 188.92 14.39 7.92 286.75 85.53 5.44 176.37 354.14 12.58 217.52 348.93 307.01 1.31 282.61 153.45 125.30 293.81 320.42 314.46 120.92 346.18 338.94 166.38 162.71 358.59 232.45 201.45 6.09 32.28 9.79 299.97 30.50 121.54 345.10 231.28 301.42 299.48 186.10 227.33 144.22 249.58 39.78 233.83 318.41

log σ −4.52 −5.00 — — — −5.47 −6.06 −4.66 −5.67 — — — −6.33 — −7.24 −5.72 — — — −3.75 −6.52 −6.49 −5.51 — — −6.10 −2.98 −4.68 — — −4.80 −4.29 −3.99 −6.45 −4.77 −5.98 −5.12 −5.45 — −5.60 −4.94 — −5.00 — −6.99 — −5.39 −6.38 −2.89 −4.36 −4.36 — −5.81 −5.45 −0.51 — — −6.03 −5.52 — −6.00

( p 1 : p2 )

Figure

2 : −3 — 4 : −7 2 : −3 2 : −3 2 : −3 — — 4 : −7 — 2 : −3 3 : −5 2 : −3 — 3 : −4 — — 2 : −3 2 : −3 2 : −3 — — 2 : −3 — — 2 : −3 2 : −3 — — — 2 : −3 2 : −3 — 6 : −11 — — 4 : −7 — — — 2 : −3 — — 1 : −2 — 2 : −3 — 2 : −3 — — — — — 2 : −3 — — — 2 : −3 2 : −3 — — 2 : −3 —

10u 11t 11p 10i 10s 10h 12k 11l 11r 11o 10b 12h 10r 11e 13a 12m 11s 10n 10p 10e 12c 12n 10j 12j 12e 10g 10o 11a 13f 12f 10t 10c 12g 12i 12o 11b 11g 11j 11n 11f 10k 12b 11h 13d 11c 10f 11k 10l 13c 11i 13b 12d 11u 10q 12a 11q 13e 10d 10m 11d 12l 10a 11m

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

21

FIG. 10. Frequency map of Kuiper Belt objects in the (a0 , ω0 ) plane. The initial values of i 0 , M0 , and Ä0 depends on the body (Table V). The color code is for log σ (units of Ma−1 ).

22

ROBUTEL AND LASKAR

FIG. 11. Frequency map of Kuiper Belt objects in the (a0 , ω0 ) plane. The initial values of i 0 , M0 , and Ä0 depends on the body (Table V). The color code is for log σ (units of Ma−1 ).

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

23

FIG. 12. Frequency map of Kuiper Belt objects in the (a0 , ω0 ) plane. The initial values of i 0 , M0 , and Ä0 depends on the body (Table V). The color code is for log σ (units of Ma−1 ).

24

ROBUTEL AND LASKAR

formed. We will give no details here for the others resonances, but the same reasoning still holds. In this region, only a few bodies are trapped in mean motion resonances with Neptune. Three of them appear to be in resonance (4 : −7). For 1997CT29 (Fig. 11g) it is difficult to say without further study whether it is in the libration island while 1994GV9 (Fig. 11r) is close to the hyperbolic orbit. 1993FW (Fig. 11p) seems to be in the small island of this resonance, while it is very clear that 1994JS (Fig. 12h) is almost in the middle of one of the three islands of the (3 : −5) resonance with Neptune. 6.3. Kuiper objects with 44 AU ≤ a ≤ 48 AU The objects with initial semi-major axis in the interval [44, 48] AU are gathered in Fig. 12i to Fig. 12o. The step size for a0 is still equal to 0.025 AU and the initial eccentricity still increases from first picture to the last. The (5 : −9) resonance with Neptune is the same in this sample of bodies as in the previous one. The other main resonances encountered are (6 : −11) at 45.1 AU, (7 : −13) at 45.5 AU, (8 : −15) close to 45.8, (9 : −17) at 46.05 AU, (10 : −19) close to 46.2 AU, and especially (1 : −2) on the right side of the pictures. Their width increases with eccentricity and covers an important area in Fig. 12o. It is interesting to note the enlargement and the modification of the shape of the resonance (1 : −2) from Fig. 12i to Fig. 12o. For moderate inclination and eccentricity (Fig. 12i to Fig. 12l), its libration island has roughly the shape expected from a simple pendulum model, but for higher eccentricity, especially in Fig. 12o for 1997CR29 (e0 = 0.213, i 0 = 17.72◦ ), the deformation due to the harmonic 2(λ8 − 2λ + $ ) = 2(λ8 − Ä − 2M − ω) is very sensitive and we obtain the characteristic shape of the pendulum with two harmonics of similar amplitude (Fig. 7). In this region, only 1997CQ29 (Fig. 12i) seems to be in the hyperbolic region of the (6 : −11) mean motion resonance. Nevertheless, it should be noted that in most of these plots, the trajectories are quite regular, and to be close to an hyperbolic orbit will not mean that the behavior of the body will be very chaotic. In particular, its wanderings will most probably still be very limited. 6.4. The Six Other Objects In Fig. 13 are gathered three objects with semi-major axis lower than 41 AU that are not in the (2:−3) resonance (Figs. 13a– 13c), and the three objects with semi-major axis greater than 47 AU (Figs. 13d–13f). The first three figures show various aspects of the inner Kuiper belt (between Neptune and 41 AU). The first picture shows a very chaotic region where 1995DA2 is inside one of the three lobes of the (3 : −4) resonance with Neptune. In Figs. 13b and 13c the diffusion rate σ is globally small, but the chains of resonant islands take up a very large surface. Thus, 1998HN151 is between two resonant chains and should have a chaotic motion (its diffusion rate is equal to σ = 10−2.89 ). It is interesting to note that in these two figures, the distortion of the (2 : −3) resonance

FIG. 13. Frequency map of Kuiper Belt objects in the (a0 , ω0 ) plane. The initial values of i 0 , M0 , and Ä0 depends on the body (Table V). The color code is for log σ (units of Ma−1 ).

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

25

seems to be in an intermediate situation between 2 and 3 lobes (see explanations above). The next figures present the same region for high eccentricity. 1997SZ10 (Fig. 13d) is very close to the edge of the (1 : −2) resonance and will probably fall into the very chaotic region outside this resonance. 1998SM165 (Fig. 13e) is itself in this large chaotic area generated by the close encounters with Neptune, and has the greatest diffusion rate of the bodies studied here, that is, σ = 10−0.51 . The last Kuiper belt object, 1996TL66 (Fig. 13f) with an initial semi-major axis of 84 AU and an eccentricity of 0.58, is itself in a chaotic region where σ is close to 10−2 . In this region where the islands are wide and close to each other, the resonances have a mutual influence that distorts these islands and leads to very dissymmetric chains. This dissymmetry can be compared with the distortion discussed in Section 3.3.2 for large semi-major axis. 7. FREQUENCY MAP OF THE v-ANDROMEDAE PLANETARY SYSTEM

Finally, in this last section, we will apply FMA to the understanding of dynamics of the newly discovered (Butler et al. 1999) v-Andromedae planetary system. In the near future, many planetary systems will be discovered. It is thus important to have a method for studying their dynamics in a fast and systematic way, one of the outcomes of these studies being the search for stability conditions that could constrain the location of eventual additional planets in the observed system. As the Doppler detection method cannot determine neither the inclination of the planetary orbital plane nor the mass of the planet, only the product of its mass by the sine of the inclination of the orbital plane to the plane of the sky (M sin(i)) can be obtained. Because sin(i) cannot be determined, only a lower limit of the planetary masses can be obtained. On the other hand, eccentricities and semi-major axes are accessible, whereas a large uncertainty persists in the angles of position of the planets. Butler et al. (1999) provides three different sets of parameters for this planetary system: the Lick Observatory data, AFOE

TABLE VI Initial Conditions of the v-Andromedae Planetary System Planet Mass (Jup.) a (AU) e ω (◦ ) M (◦ ) log σ

a 0.72 0.059 0.041 197 0 −4.99

b 1.92 0.829 0.232 195 226 −5.76

c 4.11 2.495 0.362 201 88 −5.58

Note. On the first line are the assumed masses of the planets (in Jupiter mass); on the four following are their initial conditions, and the diffusion rate (measured on a time span of 4500 years).

FIG. 14. Massless body in the v-Andromedae system: Frequency map in the (a0 , e0 ) plane for M0 = ω0 = 0◦ . The color code is for log σ , measured on a time span of 4500 years.

spectrometer data at the Whipple Telescope, and a combination of the two previous samples. Numerical integrations (Lissauer 1999) show that the second data set leads to an unstable system, while the planetary behavior is stable for several millions years for the two other samples. In our experiment, we thus assumed that the system is planar and that the masses of planets are equal to their lower limit. The other initial conditions are chosen as in the Lick data set and are reproduced from Table VI. The mass of the central star is fixed to 1.3 solar mass. The planetary system with test particles is integrated over two consecutive time intervals of 9000 years. This apparently short time, which corresponds to more than 1,400,000 revolutions of the first planet, 27,000 of the second, and 5,180 revolutions of the third, is in fact long enough to allow us to obtain a full vision of the stability regions and the location of the orbital resonances

26

ROBUTEL AND LASKAR

in the v-Andromedae system (Fig. 14). As in Section 2, the initial phase angles of the test particles are fixed to M0 = ω0 = 0, and the initial eccentricity e0 is sampled for 0 to 0.95 with a resolution of 0.05. Figure 14a is the frequency map for the region between the first two planets (0.06 ≤ a0 ≤ 0.826 AU with a step size of 0.0038 AU). The two small vertical gaps at 0.094 and 0.015 AU mark the presence of the (1 : −2) and (1 : −4) resonances with the inner planet. The three other gaps at respectively 0.283, 0.328, and 0.4 AU correspond to the (5 : −1), (4 : −1), and (3 : −1) mean motion resonances with the second planet. The (5 : −2) resonance with the second planet marks the beginning of the region where almost all test particles are rapidly ejected after close encounters with planets. In the region between the first two planets, we can see in Fig. 14a that there is only a small spot (a0 ≈ 0.25 AU, e0 ≤ 0.15) of apparently regular initial conditions where particles could survive for a long time, although for very long times, one should also check for secular resonances and chaos in this area. Figure 14b (a0 ∈ [0.83, 2.5] AU with a step size of 0.00835 AU) shows clearly that no particle can survive there, and we are in presence of a very chaotic behavior due to resonance overlap and close encounters with the planets. The situation is the same in Fig. 14c (a0 ∈ [2.5, 5] AU with a step size of 0.0125 AU). As we move away from the third planet, we find more stables areas (Figs. 14d and 14e). The region between 5 and 7.5 AU is still dominated by resonances (see Table VII) and there should not be stable regions apart from the resonant islands. Beyond 7.5 AU, the diffusion rate decreases, and we have a stable zone similar to the zone observed beyond Neptune in our Solar System. This should thus be the place to look for additional bodies. The present study was done with a massless particle, but it should be noted that for a massive planet, the chaotic regions will even be TABLE VII Mean Motion Resonances with the v-Andromedae’s Third Planet Clearly Identified in Fig. 14 a0 (AU)

p1

p2

w (AU)

4.67 5.28 5.83 6.42 6.88 7.35 7.92 8.35 9.03 9.25 9.75 10.12 10.60

2 1 2 1 2 1 2 1 2 1 2 1 2

−5 −3 −7 −4 −9 −5 −11 −6 −13 −7 −15 −8 −17

0.06 0.11 0.21 0.34 0.14 0.26 0.21 0.39 0.26 0.34 0.24 0.41 0.31

Note. a0 corresponds to the location of the base of the libration island. The integers p j define the resonance, and w is its approximate maximal width measured along a horizontal line.

larger, and the regular regions will be smaller than those ones found here. 8. CONCLUSIONS

The frequency map analysis allows us to set up rigorous tools for studying numerically in an extensive way the dynamics of particles in the Solar System. A very important aspect of this method, which is not always perceived, is the understanding that one can fix all the angle-like variables (here the initial values of the phase angles M0 , $0 , Ä0 ) to given values and still have a complete view of the dynamics of the problem. The only variations that occur for different choices of the phase angles are differences in the crossing of the resonant regions (Section 4). The other fundamental feature is the possibility of obtaining a diffusion indicator in a very short integration time. Here σ is a measure of the diffusion of the proper mean motion, which we related to a diffusion in semi-major axes. This very short time of integration allowed us to span a huge numbers of initial conditions and thus to obtain a complete view of the dynamics of the system. In each of the figures 2 and 9 are summarized, for example, the results of the numerical integration of 192,000 test particles, and for the whole paper, we integrated 2,048,000 orbits. Of course, this was only possible because we used a very short time of integration about 1000 time smaller than what is often used now for similar studies. The next feature of FMA is the possibility to detect in a systematic way all the resonances that are of some importance in the dynamics. Therefore, not only do we know the chaotic areas and the regular regions, but we also have all the resonant islands with their true shape. We can thus also understand the development of chaotic regions arising from resonance overlap. We have obtained here a complete dynamical map of the Solar System (Fig. 2, Fig. 9), which can be used as a framework for more specific issues, for which a more detailed analysis could be done. This gives in particular a striking view of the outer planetary system (the Kuiper belt), which is clearly separated in four components: the collisional region, the resonance overlap region, and the stable region, linked together by the resonant islands (Section 3.4). FMA also allowed us to analyze in detail the geometry of the resonance in the complete phase space (Section 4), and we showed that a simple pendulum model of resonance can predict very well the location of the islands, while the shape of the island sometimes requires adding a few harmonics. We could thus give for all Kuiper belt objects with sufficiently well-known orbital elements an identification picture of the global dynamics in their vicinity. This is particularly important as the orbital elements of these objects are not very well determined, so the knowledge of a single orbit is not sufficient for understanding the dynamics of the true object. Another important aspect of the present study is the setting of an effective tool aimed at analyzing the numerous planetary systems that will certainly be discovered in the near future. As

FREQUENCY MAP AND GLOBAL DYNAMICS IN THE SOLAR SYSTEM I

an illustration, we presented the analysis of the dynamics of the newly discovered v-Andromedae system, which shows the location where additional planets could exist.

27

Dumas, S., and J. Laskar 1993. Global dynamics and long-time stability in Hamiltonian systems via numerical frequency analysis. Phys. Rev. Lett. 70, 2975–2979. Duncan, M., H. Levison, and S. Budd 1995. The dynamical structure of the Kuiper belt. Astron. J. 110, 3073–3080.

8.1. Perspectives Some issues still need to be discussed. For orbits with very small diffusion (σ < 10−6 , i.e., the dark blue region in Fig. 2), we are probably ensured of the stability over very long time, similar to the age of the Solar System. On the other hand, orbits with σ > 10−3 (from red to yellow–green) are chaotic with macroscopic instabilities, and are for particles that will be expelled during the lifetime of the Solar System, but uncertainty remains for the particles in the intermediate situation (10−6 < σ < 10−3 ). In this case further study is necessary, as different options can take place: (i) This value is overestimated due to uncertainty of the FMA method. This is particularly true in the vicinity of the resonances, where there is a deformation of the frequency map. This case can be ruled out by a longer integration as the FMA accuracy increases as 1/T 4 (see Eq. (10)). (ii) The value of σ is properly determined and the chaotic behavior comes from mean motion resonances. In this case, the diffusion of the particle can be understood by the present study. (iii) The value of σ is properly determined, but this diffusion results in fact from chaotic behavior in the secular system. In this case, a further study is needed as the diffusion in the secular system can be far more important than in the short periods. This could lead to large increase of eccentricities, and thus very unstable motion. This will also require a longer integration time to be fully understood, as we will search in this case to determine properly the secular frequencies. This should happen for example in the zones under the collisions lines in between Mercury, Venus, the Earth, and Mars. We are thus still left with uncertainties, but these are limited to specific areas of the full map, where more detailed analysis could be conducted, after having thus restricted very much the set of parameters of the system. ACKNOWLEDGMENTS The computations have been made on an IBM-SP2 at CNUSC/Cines. A contribution from PNP-CNRS is acknowledged.

REFERENCES Arnold, V. 1988. Dynamical Systems III. Springer-Verlag, Berlin. Beaug´e, C. 1994. Asymmetric libration in the exterior resonances. Celest. Mech. Dynam. Astron 60, 255–248. Butler, P., G. Marcy, D. Fischer, T. Brown, A. Contos, S. Korzennik, P. Nisenson, and R. Noyes 1999. Evidence for multiple companions to upsilon Andromedae. Astrophys. J. 526, 916–927. De la Barre, C., W. Kaula, and F. Varadi 1996. A study of orbits near Saturn’s triangular Lagrangian points. Icarus 121, 88–113.

Evans, W., and S. Tabachnik 1999. Possible long-lived asteroid belts in the inner Solar System. Nature 399, 41–43. Ferraz-Mello, S. 1994. Dynamics of the asteroidal 2/1 resonance. Astron. J. 108, 2330–2337. Gladman, B., and M. Duncan 1990. On the fates of minor bodies in the outer Solar System. Astron. J. 100, 1680–1693. G´omez, G., A. Jorba, J. Masdemont, and C. Sim´o 1993. Study of the Poincar´e maps for orbits near Lagrangian points. ESA Technical report. Grazier, K., W. Newman, F. Varadi, W. Kaula, and M. Hyman 1999a. Dynamical evolution of planetesimals in the oute Solar System I. The Jupiter/Saturn zone. Icarus 140, 341–352. Grazier, K., W. Newman, F. Varadi, W. Kaula, and M. Hyman 1999b. Dynamical evolution of planetesimals in the oute Solar System II. The Saturn/Uranus and Uranus/Neptune zones. Icarus 140, 353–368. Henrard, J., and A. Lemaitre 1983. A second fundamental model of resonance. Celest. Mech. Dynam. Astron. 30, 197–218. Holman, M. 1997. A possible long-lived belt of objects between Uranus and Neptune. Nature 387, 785–788. Holman, M., and J. Wisdom 1993. Dynamical stability in the outer Solar System and the delivery of short period comets. Astron. J. 105, 1987–1999. Laskar, J. 1990. The chaotic motion of the Solar System. A numerical estimate of the size of the chaotic zones. Icarus 88, 266–291. Laskar, J. 1993. Frequency analysis for multi-dimensional systems. Global dynamics and diffusion. Physica D 67, 257–281. Laskar, J. 1997. Large scale chaos and the spacing of the inner planets. Astron. Astrophys. 317, 75–78. Laskar, J. 1999. Introduction to frequency map analysis. In NATO ASI Hamiltonian Systems with Three or More Degrees of Freedom (C. Sim`o, Ed.), pp. 134– 150. Kluwer, Dordrecht. Laskar, J. 2000. On the spacing of planetary systems. Phys. Rev. Lett. 84, 3240– 3243. Laskar, J., and D. Robin 1996. Application of frequency map analysis to the ALS. Particle Accelerator 54, 183–192. Laskar, J., and P. Robutel 1993. The chaotic obliquity of the planets. Nature 361, 608–612. Laskar, J., and P. Robutel 1995. Stability of the planetary three-body problem. I Expansion of the planetary Hamiltonian. Celest. Mech. Dynam. Astron. 62, 193–217. Laskar, J., C. Froeschl´e, and A. Celletti 1992. The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping. Physica D 56, 253–269. Levison, H., and M. Duncan 1993. The gravitational sculpting of the Kuiper belt. Astrophys. J. Lett. 406, 35–38. Levison, H., and M. Duncan 1994. The long-term dynamical behavior of short period comets. Icarus 108, 18–36. Levison, H., and M. Duncan 1997. From the Kuiper belt to Jupiter-family comets: The space distribution of ecliptic comets. Icarus 127, 13–32. Lissauer, J. 1999. Three planets for Upsilon Andromedae. Nature 398, 659–660. Luu, J. 1994. The Kuiper belt, In symposium IAU 160 (A. Milani, M. Di Martino, and A. Cellino, Eds.), pp. 31–44. Kluwer, Dordrecht. Malhotra, R. 1996. The phase space structure near Neptune resonances in the Kuiper belt. Astron. J. 111, 504–516. Message, P. 1958. The search for asymmetric periodic orbits in the restricted problem of three bodies. Astron. J. 63, 443–448.

28

ROBUTEL AND LASKAR

Milani, A., and Z. Knezevic 1990. Secular perturbation theory and computation of asteroid proper elements. Celest. Mech. Dynam. Astron. 49, 347–411. Moons, M., and A. Morbidelli 1993. The main mean motion commensurabilities in the planar circular and elliptic problem. Celest. Mech. Dynam. Astron. 57, 99–108. Moons, M., and A. Morbidelli 1995. Secular resonances in mean motion commesurabilities: The 4/1, 3/1, 5/2 and 7/3 cases. Icarus 114, 33–50. Morbidelli, A. 1997. Chaotic diffusion and the origin of comets from the 2/3 resonances in the Kuiper belt. Icarus 127, 1–12. Morbidelli, A., F. Thomas, and M. Moons 1995. The resonant structure of the Kuiper belt and the dynamics of the first five trans-neptunian objects. Icarus 118, 322–340. Nesvorny, D., and S. Ferraz-Mello 1997. On the asteroidal population of the first-order jovian resonances. Icarus 130, 247–258. Nesvorny, D., and A. Morbidelli 1998. Three-body mean motion resonances and the chaotic structure of the asteroid belt. Astron. J. 116, 3029–3037. Papaphilippou, Y., and J. Laskar 1996. Frequency map analysis and global dynamics in a two degrees of freedom galactic potential. Astron. Astrophys. 307, 427–449. Papaphilippou, Y., and J. Laskar 1998. Global dynamics of triaxial galactic models through frequency map analysis. Astron. Astrophys. 329, 451–481.

Poincar´e, H. 1902. Sur les plan`etes du type d’H´ecube. Bull. Astron. 19, 289–315. P¨oschel, J. 1982. Integrability of Hamiltonian systems on Cantor sets. Comm. Pure Appl. Math. 25, 653–695. Robutel, P. 1993. Contribution a` l’´etude de la stabilit´e du probl`eme plan´etaire des trois corps. Th`ese de l’Obsevatoire de Paris. Sim´o, C., G. G´omez, A. Jorba, and J. Masdemont 1995. The bicircular model near the triangular libration points in the RTBP. In From Newton to Chaos: Modern Techniques for Understanding and Coping with Chaos in n-Body Dynamical Systems (A. E. Roy and B. A. Steves, Eds.), pp. 343–370. Plenum, New York. Standish, E. M. 1990. The observation basis for JPL’s DE200, the planetary ephemerides of the Astronomical Almanac. Astron. Astophys. 233, 252– 271. Torbett, M. 1989. Chaotic motion in a comet disk beyond Neptune: The delivery of short-period comets. Astron. J. 98, 1477–1481. Wisdom, J. 1980. The resonance overlap criterion and the onset of stochastic behavior in the restricted three-body problem. Astron. J. 85, 1122–1133. Wisdom, J. 1983. Chaotic behavior and the origin of the 3/1 Kirkwood gap. Icarus 54, 51–74. Wisdom, J., and M. Holman 1991. Symplectic maps for the n-body problem. Astron. J. 102, 1528–1538.