Relaxation-chaos phenomena in celestial mechanics

Physica 26D (1987) 85-122 North-Holland, Amsterdam

RELAXATION-CHAOS PHENOMENA IN CELESTIAL MECHANICS I. ON WISDOM'S MODEL FOR THE 3 / 1 KIRKWOOD GAP

Jair KOILLER*, lnstituto de Matematica, UFRJ, Caixa Postal 68530, Rio de Janeiro, R J, Brazil, 21944

J.M. BALTHAZAR and T. YOKOYAMA Departamento de Matematica Aplicada, UNESP, Caixa Postal 178, Rio Claro, SP, Brazil, 13500 Received 18 June 1986 Revised manuscript received 28 October 1986

Sudden eccentricity increases of asteroidal motion in 3/1 resonance with Jupiter were discovered and explained by J. Wisdom through the occurrence of jumps in the action corresponding to the critical angle (resonant combination of the mean motions). We pursue some aspects of this mechanism, which could be termed relaxation-chaos: that is, an unconventional form of homoclinic behavior arising in perturbed integrable Hamiltonian systems for which the KAM theorem hypothesis do not hold.

1. Introduction

Since Kirkwood's 1867 discovery [1] of gaps in the distribution of the semimajor axis of the asteroids in the belt between Mars and Jupiter, explaining this phenomenon has been considered a major challenge in Celestial Mechanics. A review of different theories presented up to 1986 is given in [2]. Recent numerical and analytical work by J. Wisdom [3] indicates that purely gravitational effects may account for the gaps in several of the resonances. First, using an ingeniously devised algebraic mapping instead of the standard ODEs solvers, he was able to integrate the (high-frequency terms removed) equations near the 3/1 resonance, for a much larger time span than previous works. Surprisingly enough, sudden eccentricity jumps from - 0 . 1 to - 0 . 3 5 appeared for typical trajectories. This discovery introduced "a new possibility for the origin of the Kirkwood gaps. The hypothesis is that asteroids near the commensurabilities undergo larger jumps in eccentricity, thus becoming Mars crossers. They are subsequently removed from the gaps through collisions or close encounters with Mars" [3a]. These eccentricity jumps were confirmed later by Murray and Fox [4], integrating full and averaged equations of motion by reliable (although more expensive) ODEs solvers. Recently, Wisdom has also produced a theoretical explanation for the phenomenon [3c]. It is based on the observation that the critical angle ~p (three times the mean longitude of Jupiter minus the mean longitude of the asteroid) and its conjugate momentum • (essentially a function of the semimajor axis a) oscillate in a faster rate than the eccentricity e and longitude of periapse w. The motion of (cp, ~) can be mathematically identified to that of a pendulum with slowly varying coefficients (functions of e, w). Assuming the "adiabatic invariance" of the action I of this pendulum, averaged equations for (e, w) are obtained. However, the procedure breaks down when the rates of (rp, ~b) and (e, w) become compatible, i.e., when the pendulum motion approaches *On a visit to the Department of Mathematics, Yale University, 1985/1986, under a CAPES/Brazil fellowship.

0167-2789/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

86

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

the separatrix. Wisdom shows that precisely this break down is the cause of the eccentricity increases: a trajectory stays in the ~0-1ibration or circulation regime with constant I until it reaches the "uncertainty zone". After some unpredictable motions inside this region, the trajectory eventually leaves it in one of the two regimes with a different action. The aim of this paper is to pursue some theoretical aspects related to [3c]. In a sequel, we will concentrate mostly on numerical experimentations. Our main results here are the following: i) We formalize the canonical transformation method so that the first perturbative term to the averaged system can be easily calculated. ii) Using the linear oscillator approximation for the libration regime we give a rough estimate of the maximum eccentricity attained and the amplitude of ~-motion (the procedure is simple and we plan to apply it to other resonances). iii) We show that there are eccentricity increases which take place within the regime of very small librations (thus occurring without a correlated jump in the action) due to a singular form of Poincarr's homoclinic bifurcation [5]. Although this is a secondary source of chaos, it arises around Hill's periodic orbit [6] and therefore may have at least some historical interest. Here eccentricity jumps to about 0.24. iv) Using Melnikov's method. Wisdom's criterion for chaos is recovered and interpreted in terms of transversal intersection of stable and unstable manifolds. v) We show that the regular homoclinic phenomenon is present when the eccentricity of the smaller primary is taken as the small parameter. In addition, we present some speculative arguments for the asymptotic behavior of the uncertainty zone width and a tentative numerical scheme to simulate uncertainty zone motions, using "Chirikov's trick" [7] in the same line of ideas of Wisd0m's mappings. Hopefully we will be able to follow this program for the sequel paper. Finally, we wish to explain why we propose the name "relaxation-chaos" for Wisdom's mechanism. Maximum eccentricities correspond to the smallest possible libration amplitude I * compatible with the constraint that (e, w) motion reaches the uncertainty zone (the circulation regime produces motions in the region inside the UZ central curve). The alternation of paths with endpoints at the UZ vaguely resembles van der Pol's famous "relaxation-oscillations" [8]. Although the analogy cannot be taken literally (since here the singular parameter/~, mass ratio of Jupiter and the Sun, appears multiplying a second derivative), it was the source of our preliminary insights before learning Wisdom's explanation. An announcement of this work appeared in [9].

2. The constrained systems 2.1. Wisdom ' s equations for the 3/1 resonance Following [3], consider the two degrees of freedom Hamiltonian system H --- ½a¢ 2 +/~R,

(1)

R = F ( x 2 + y2) + f i e j x - [ C ( x 2 - y 2) + Dejx + EejZ)] cos ~ - [2Cxy + Dyejlsincp. Here ¢p is the critical angle (three times the mean longitude of Jupiter minus the mean longitude of the asteroid); • its conjugate variable (function of the asteroid semimajor axis a, zero at exact resonance);

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

87

(x, y ) are Poincar6 variables x = ( 2 p ) 1/2 cos w,

y = (2p) 1/2 sin w,

p = ( i x l a ) l / 2 [ l -- ( 1 - - e2) 1/2] = ( I x l a ) 1 / 2 e 2 / 2

(/.tl= 1--/~),

(2)

where w is minus the longitude of the asteroid periapse, e its eccentricity. The time scale t is such that Jupiter's period is 2~r. More precisely, units are such that aj=l,

G(M+mj)=I,

nj=[G(M+mj)/a3j]l/2=l.

The parameters are given by a = - 12.98851, C = 0.8631579,

F = -0.2050694, i f = 0.1987054, D = -2.656407, E = 0.3629536.

(3)

One can check that exact resonance occurs for a * = 0.4805968 [4]. In the sequel we will leave free the parameters ej (Jupiter's eccentricity) and"/~ = m j / ( M + m j) (mass ratio of the primaries), but substitute, when convenient ej = 0.048,

/~ = 1/1047.355,

#x = 1 - / ~ .

(4)

Remark. Deriving (1), a deceivably simple looking equation which couples a pendulum to an harmonic oscillator, was already a major achievement. For details and more precise astronomical definitions, see [3a].

2.2.

A singular perturbation problem

Performing the time-scale change ~"=/~t the equations of motion become ( ' =

xt=-Ry,

d/d,r)

y'=Rx,

(5a)

#q~" = - a R ~ .

(5b)

Notice that r is the natural scale for (e, w) motion. Setting ~t = 0 we obtain a since (5b) becomes

= 0.

"constrained system",

(5c)

We will see on section 3 that (5a, 5c) approximate only very particular classes of motions of (1). Nevertheless this study is almost effortless, and our view is that it should be done as a starting point, as we now proceed to do for the 3 / 1 commensurability, R ~ = A s i n cp -

B cos ¢p,

A = C(x 2 _ y 2 ) + Dejx + Eel,

B = 2Cxy

+ Oyej.

(6)

88

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

Solving for ¢p = ¢p(x, y) we obtain, provided A 2 + B 2 :g 0. branch 1:

cos ~o(1) = A l P ,

branch 2:

cos ¢p(2)= _ A L P ,

sin (]9(1)

=

B/P,

p2 = A 2 + B 2,

(7)

sin ¢p(2)= _ B / P .

Inserting (7) into (5a) we obtain the constrained systems in the (x, y) plane (upper signs refer to branch 1):

x' = - 2Fy -T-2CyA/P + (2Cx + D e j ) B / e ,

(8a)

y ' = 2Fx + ffej • (2Cx + D e j ) A / P T- 2CyB/P.

(8b)

These equations are of Hamiltonian type, with Hamiltonian functions R(X),(2) = F ( x 2 + y 2 ) + ffejx • P.

(9)

This follows from the fact that if/~(x, y) = R(x, y, ~(x, y)) then

k x=R x+Rr%=Rx,

etc.

whenever cp(x, y) is a solution of (5c).

Remarks. i) Fixing (x, y), (5b) becomes a one degree of freedom system. Then ¢p(1) is an unstable while stable equilibrium point. Mathematically, (5b) is identified with the equations for a pendulum ~ p " = a/t- - P s i n ~ ,

tp=q0-~(E)(x,y)

(a<0).

q9 (2)

is a

(10)

ii) In common applications of constrained, or implicitly defined ODEs (see [10abc]), the singular parameter appears on a first derivative and the constraint defines an attractor set. Offhand, one can only hope that (5a, 5c) approximate some of the motions (we thank P. Holmes for stressing this point); see, however [10d] for a new approach that promises some surprises.

2.3. The constrained systems Level curves of P(x, y) are depicted in fig. 1. Notice that P = 0 at

Xl' x2 =

- D +_~/D 2 - 4CE 2C e j.

(11)

The equilibrium points at the x axis are easily identified: (8a) is automatically satisfied and (8b) yields

2Fx + f f e j + (2Cx + Dej) sign A

0.

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

A>O B
A
~ %

%@

A>O B>O

A O ,, I s

/ A O

X A
89

A>O B>O

~>o B
Fig. 1, Level c u r v e s of P (sketch).

Solving for x we get

+_D-ff x ! - 2( F_T_C) er

(12)

W e also designate by D X0 =

(13)

-- ~-~ej

the s y m m e t r y center of the equilateral hyperbola A = 0 and the intersection of the two lines B = 0. Substituting the numerical values we obtain x I = 0.1433eJ,

x 2 = 2.9342eJ,

x + = 1.3364eJ,

x

x o = 1.5387ej,

=1.8625eJ.

(14)

Since x I < x+ < x 0 < x _ < x2, and since xz, x 2 are the roots of A ( x , 0 ) = 0, a quadratic equation, we have sign A < 0 at x +. Hence:

Proposition 1. Branch I has the equilibrium point x s = x_, y = 0 while b r a n c h 2 has the equilibrium point x n = x+, y = 0. Both are unstable.

Proof. It remains to verify only the last assertion. A tedious direct calculation gives R ~ ) ( X s , 0 ) = 1.3, R~)(XH,0)=--2.1,

RCy~(Xs,0) = - 2 . 2 , R(y~(XH,0)=I.5,

R~y)(Xs,0) = 0, R~y)(XH,0)=0.

(15) []

Remark. The critical points (x +, 0) yield equilibrium solutions of the unconstrained system (5ab) as well

90

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

I

j, a,"

s/

J

•

s

1 • %

*

~

,,, x:

x:

--"~x ITx [2:

/Ill//ill

R(2)

N

%%

/-

s~

\

Fig. 2. Graphs of R(1)(x,0) and R(2)(x,0) (sketch). For the energy level shown, only a small region around x 2 is prohibited a priori (see section 4).

precisely those known in the hterature as Hill's [6] and Sinclair's [11]. From (2),

ell=e+

_

x+ (~1a,)1/4

=0.077,

e s= e

-

X

(~1a.)1/4

=0.107.

(16)

Notice that (we now change notation, H = Hill, S = Sinclair)

~(2)(x.,O) = 0,

~")(Xs,0) =,,.

Proposition 2. The solution curves emanating from (x s,0) for R (*) and (x H, 0) for R (2) form "figure eight" loops.

Proof. This follows immediately from plots of R(1)'(E)(x,0) versus x (Fig. 2). These can be sketched by elementary algebra. In fact,

R(a)'(2)(x,0) = Fx 2 + ffejx -T-[Cx2 + Dejx + Ee][, or, equivalently, (F-T-C)x

R(')'(2)(x'0)=

2 +(F-T-D)e,x-T-Ee 2, x ~ ( x , , x 2 ) ,

( F + C ) x 2 + ( f f ++_D)eJx +_Ee],

x~(x,,x2),

(17)

where we recall that upper signs refer to branch 1. Thus the extreme x values for the loops are the roots of

91

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

the quadratic equations

( F - C)x 2+ ( ff - D ) e s x - Ee 2= R~l)(xs) = [ - ( f f + D ) 2 / a ( F +

C) + E]e 2,

(F+ C)x 2 + (P+D)ejx + Ee~ = R(2)(xH) = [ - ( P - - O ) 2 / 4 ( F - C) - E] ej2. We obtain (cf. fig. 2) x~l)=-0.a7ej,

x(41)=3.14ej,

x3(2)=-0.a3ej,

(18a)

X(n2)=n.17ej,

corresponding to the eccentricity values (we use e = [xl/(tqa *)1/4) e~1) = 0.027,

e4(1)= 0.181,

e~2) = 0.025,

(18b)

e4(z)= 0.241

[] We also compute, for later use: R(1)( x 1, 0) = R (2)(x 1, 0) = 0.02e 2, R¢I)(x2,0) = R(2)(x2,0) = - 1.18e 2, (19)

R~l)(Xs,0) -- - RO)ix¢l) - 1.93e 2, ~, 3 , 4 ,0 ]~ = R(2)(XH,0) -R(2)[x(2),0 ]] = 1.54e 2. -~, 3 , 4

Proposition 3. Far enough outside the loops, the contour lines are approximately circular, centered at (Xn,0) for R ~1) and (Xs,0) for R ~2) (notice the reversal of x n and Xs).

Proof. We first expand p 2 = h 2 + n 2 = C 2 ( x 2 + y 2 ) + 2 C D e j x ( x 2 + y2) + e][ D 2 ( x 2 + y2) +

2EC(x 2 - y2)]

+ 2DExe 3 + E2e 4.

(20)

Equivalently P = C(x 2 +y2)(1 + u) 1/2 where ej

U-C(xE+y2)

(

[D22E(x2-y2)] 2Dx+ej --C+ C(xE+y2)

2DEx E2e 3 +C(x2+y2)e~+c(x2+y2)

)

(21)

.

For lul << 1 we may approximate (1 + u) 1/2 - 1 + u/2 so

p _ C(x 2 + y2) + Dejx,

(22)

where we have dropped out the O(ej2) terms. Substituting into (10) we get R (1)'(2) = ( F ~

C)(x 2 +y2) + (ff+D)ejx

+ 0(e2),

(23)

provided (this condition comes from the requirement that lul < 1). ( X -t- D e j / C ) 2 _[_y2 >

D 2e2/C.

(24) []

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

92

3. Exponentially small homoelinic phenomena 3.1. Transformation to action-angle variables. The averaged system We now present a formalization of the canonical perturbation method used by Wisdom [3c, appendix]. Consider the auxiliary one degree of freedom system for variables (% A)

h ° = ½aA2 + R(x, y, ~p),

(25)

where x, y are parameters and /, has been dropped out. Pick a parametrized family of action-angle coordinates defining a regime of ¢p-motion (say, librations or circulations)

¢p=cp(I,O;x,y),

A=A(I,O;x,y),

h°=h°(I,x,y).

It is well known that this canonical transformation can be obtained from the generating function

a ( I, cp; x, y ) = via 0 =

GI,

(2,(x,y)t(2/a)[h 0 ( I , x , y ) - R(x, y,~p')] dcp'

(26)

A = G~ and where h°(I, x, y) is the inverse function of

I = I ( h ° ; x , y ) = ~ - ~l ~ c ~ V / ( 2 / a ) [ h ° - R ( x , y , ~ ) ] d%

(27)

o

Here Cho is a family of closed level curves of the "rescaled" Hamiltonian (25). Allowing now (x, y) to recover the status of canonical conjugate variables, consider the full generating function

F= x~ + eG( I, ~, x, ~),

(28a)

for the transformation of variables (if, f, I, 0) ~ (x, y, ~, ~) given by

y=rx,

~=Fy,

e o = r t,

Cb=F~,

or

~=x+~6~, y=y+~Gx, o=a,, ~=~G~. If we choose e = ~

(28b)

then the Hamiltonian becomes

/-/= ~ [ ~ . ( * / ~ ) 2 + R(x, y , ~ ) ]

= ~ [½.G~ + R(x, y + ~G~, ~)] = . [(~.G~ + R(x, ~, ~)) + ~R.Gx + ~ ( ~ ) ] = ~[h°(I, x, y) + ~R,Gx + 0(~:)]

=/z [h°(I, Y, .~) + ehl(ff, fi, I, 0) + O(e2)],

(29a)

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

93

where

h 1 = RyG x - h°Gy

(29b)

can be calculated at (£, Y, I, 0) since the error may be transported to the d~(e2) term. Here cp is replaced by cp=~p(Y, Y, 1,0) inverting 0 = Gi(I,~p, ~, 3). (Actually any number of terms of the series H = /~[E, > 0h"(~, Y, I, 0)] can be computed, but the expressions become increasingly complicated). It should be stressed that the symplectic form [12a] now acquires a singular form dxdy+

d~dep=dYdfi+edldO,

(30)

so the equations of motion for I, O have a factor 1/e. In the intermediate time scale o = et it follows that d Y / d o = - ~0 [ehO(i,y~,3)+O(e2)] d 3 / d o = O [ehO(i,~,y)+O(e2)] ' (31)

dI/do = - h[eht(2'

Y, 1 , 0 ) + O(e2)],

dO~do = O-~ [ h ° ( I ' 2, .~) + d~(e)]. In a domain such that ~0(I, x, y) = (Oh°/OI)(I, x, y) is bounded away from zero, the averaging theorem [12a, $52] can be applied:

Theorem 1. Solutions of (31) differ by less than ce in a time interval 0 < o < T/e from solutions of the averaged system eh°(I, x, y). In particular, I is an adiabatic invariant. Remarks. i) The constant c depends on bounds on C 2 norms of h ° and h 1, but it is independent of e. It is directly proportional to T and it is inversely proportional to the lower bound on I,ol. See the proof in [12, $52]. ii) As the "slowly varying parameters" x, y are intrinsically coupled with (% ~), the adiabatic invariance of ! seems at first sight contradictory to Arnold's phrasing of "physicist's definition": "the person changing the parameters of the system must not see what state the system is in" [12]. Actually, the action invariance follows precisely from the Hamiltonian character of (31): namely, the average of 3h/O0 is zero. Example. The regime of small librations. Using the harmonic oscillator approximations, elliptic functions can be avoided and we have R = R(2)(x, y) + ½aA2 - (P/2)(cp - ¢p(2)(x, y ) ) ,

(32a)

~p = ¢p(2)(x, y) + [a/ P lX/4 22~sin O,

(32b)

a = IP/all/42¢~cos O,

(32c)

h°(I, x, y ) = R<2)(x, y) -laPll/2I,

(32d)

2 x/2f~ ~lotpll/Zi_(P/2)(ep_~p(V)(x ' y))2 dcp. G = la l "%<2)(x,y)

(32e)

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

94

Evaluating this integral, we obtain G=G(x,y,Z,cp)=$(8+sinBcos8), e = arcsin

lw41’4(‘p - 4*)(x, J4

m

(33)

.

Remarks.

- i) W. Sessin called our attention to the fact that the transformation of variables (x, y) * (x, y) is singular for Z = 0. Although the averaged system is defined for Z = 0, h”(O, x, y) = R(*)(x, y), eq. (33) shows that the perturbative term h’ blows up like l/a. ii) In order to study the behavior of the transformation of variables near the separatrix, observe that

Gx=

$$ $!+-dho(l, x,y)-R(vw’)l G’- W~)[~0-R(w~~2)1

(G, analogous).

(x,

The integrals will diverge only for ‘p = @(x, y). Following Wisdom [3c], a correspondence y) is set fixing B = 0, i.e., cp= (p(*)(x, y), so that

(x, Y) *

In, particular, at the separatrix level, ho = R(x, 7, C&X, j)) and hence

(284 3.2.

Exponentially small ergodic zones

We assume as a working hypothesis that for Z * f 0 small enough the unaveraged dynamics Z(t) does not leave a small neighborhood of Z* for a longer time span than guaranteed by theorem 1. We replace (31) by the simplified system (now written in the scale r = pt) d2 -= dr dY d7= d6’ -= d7

-~[hO(z*,~,~)+eh’(P,B,z*,8)], a[,o(z*,x,y)+Eh1(X”,P,Z*,8)], a2 $(z*,

2, y”),

(34)

0 = %P/az.

Current understanding [13] is that for small E # 0 “exponentially small horseshoes” are present, so the ergodic regions around the unperturbed separatrices have width O(exp( - const/e)). The fact that the frequency depends on x, y is a minor difficulty [14]; the more serious technical point is associated to the undefined status of the relations between averaging and the “Melnikov method” (see [15]).

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

95

Theorem 2. Poincarr's homoclinic phenomenon is present around Hill's loop for 1 = 0 and its deformations for small I. Although the width of the ergodic region is asymptotically exponentially small for I 4: 0, it is enhanced as I ---, 0. "Proof". We just observe from (32d) that separatrix loops of h(I*, x, y), 1" small, can be thought as deformations of Hill's loop. The previous remark (i) suggests that among these ergodic zones, that one around Hill's loop should be the most visible one. [Further work is necessary to clarify this case, but we speculate that perhaps the width is here linear with respect to el. [] We want to stress that the existence of exponentially small ergodic zones was predicted by Wisdom. They seem to be a rather subtle phenomenon. Accordingly, "I did not find chaotic zones near the unstable equilibria of the guiding trajectories, but I expect that with sufficient magnification this would be possible" [3c, $VI]. One may ask what is the effect of the "high frequency terms" [3, 4] of Jupiter's disturbing function which were truncated in the derivation of (1). We have the following (somewhat disappointing) result:

Theorem 3. Perturbing the averaged system h°(I, x, y) with high frequency terms, ergodic zones are produced around the separatrices, but the width is 0 ( e x p ( - c o n s t / / Q . Hence the effect is negligible with respect to that given by the previous theorem. "Proof ". The reader should bare the following heuristic considerations, which we think may be suggestive. The high frequencies are multiples of Jupiter mean motion and the amplitudes are proportional to/~. We can represent the smaller harmonic perturbation (using the ,r-scale)

h = h°(I *, x, y) +p(x, y, I*, 0) cos (~'//~).

(35)

Notice that the frequency is much higher (O(1/e) times) that of 0 (cf. (34)), but the amplitude is of the same order of the Harniltonian h. How does one justify averaging p before ehX?As far as we know, there is no rigorous mathematical treatment as yet available for systems like (35) (for another example, see [16]). Nonetheless, emulating Wisdom's mappings derivation [3], we use the principle that if truncating is OK, so is "untruncating". We thus proceed as follows: Let I = 0. Break R(Z)(x, y) into two parts,

hi=Rsec=F(x2+y2)+ffejX,

hn=P(x,y )

(36)

and consider a superposition of perturbing terms of the same amplitude P cos n(T//~). Chirikov's trick [7] yields the integrable Hamiltonian with "kicks" h = h I + hti82~('r//~).

(37)

The evolution of this Hamiltonian is described by the discrete mapping corresponding to the time interval Ar = 2¢r#,

(38)

obtained by superposing to the h I flow an instantaneous kick A~-(-py, p~). More precisely, we define the

96

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

algebraic mapping Xa/2

=

--

f f e j / 2 F + cos (2FA'r)(x 0 + f f e j / 2 F ) - sin (2FAz)Y0,

(39a)

Yl/2 = sin (2FA~')(x 0 + f f e j / 2 F ) + cos (2FA'r)y0; Xl = X1/2 -- A'r P y ( X l / 2 , Y t / 2 ) ,

(39b)

Yl =Yl/2 + A ' r P x ( x t / 2 , Yl/2)" [Remark. In order to make (39b) closer to a symplectic mapping (as it should) one may need to use a better approximation to the hti flow for time At.] Therefore, (39) can be interpreted as a "clumsy" way of discretizing R, since the Poisson bracket (h I, hit } ~ 0 and the flows do not commute. A similar mapping composition was studied by Aref [17] on an hydrodynamical context. He considers the motion of a marker particle on a planar inviscid flow generated by two stirrers, turned on and off at regular intervals T. If T = 0, i.e., the stirrers are operating continuously, there is a figure eight loop in the marker's phase plane, which unfolds into an ergodic zone for small T 4=0. According to Aref, " . . . the size of the ergodic region goes to zero as e x p ( - c o n s t / T ) as T ~ 0... The onset of chaos [is] a realization of the numerical instabilities to which the leap frog scheme is succeptible at large stepsize". [] Remark. In our first insight into the 3/1 resonance problem, we erroneously believed that the source of the global chaotic behavior were bursts of rapid variation of cp, making the system flip randomly between branches 1 and 2. It turns out that this idea represents sort of a "caricature" of the correct explanation, discovered by Wisdom himself, which we describe in the next paragraph. Nonetheless, we believe that it is worth to state now a mathematical problem motivated by our (wild) first guess, which is also connected with Aref's stirrers. For large eccentricities, branch 1 and branch 2 solutions are approximately concentric families of circles (with different centers). The system is geometrically identical to Aref's advection by "pulsating" stirrers. Now, if one considers the time pulses and direction of motion as random variables (satisfying certain "rules", or statistical distributions), the resulting advection behaves like a random walk, except that the underlying lattice does not exist anymore (because the flows do not commute). One may ask, however, similar questions to those arising in potential theory, (as probabilities of escape, for example). In the same line of thought, one may consider the abstract mathematical problem of random walks following noncommuting vectorfields. In the Lie group theoretical context, the prototype examples are the three dimensional Heisenberg group and the rotation group SO(3). The random walk would be along two noncommuting left invariant vectorfields. We speculate that there should be a connection with the "potential theory" of the associated hypoelliptic operators. The authors ignore if this question belongs to a well-developed mathematical topic.

4. Relaxation chaos

4.1. Wisdom ' s explanation One must remember that both the change of variables (28) and the averaging procedure (Theorem 1) break down when (Y(o), ~(o)) approaches a point where 2~rI is the area inside the separatrix of the

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

97

0.1

0.0

-0.1

-0.2 -0.1

0.0

0.1

0.2

0.3

Fig. 3. Wisdom's diagram of averaged motions with H = - 4 . 2 2 × 10 -6 (reprinted with permission).

instantaneous (~p, A) system (25), hereafter called the "frozen pendulum" (FP). The energy of the FP separatrix can be obtained using the FP unstable critical point where A = 0 and q9 = q~tl~(~, 5). We obtain h = R(I~(~, y ) as the break down condition, putting into evidence the importance of branch 1 solutions. The crux of Wisdom's method is embodied in the

Classical Uncertainty Principle. "The [adiabatic] approximation will also break down in some neighborhood of these curves as well..., but its actual extent remains undetermined. This region will be called ' the zone of uncertainty' [UZ]. It is in this zone that the actual trajectories can no longer be followed by the averaging method. The zone of uncertainty separates regions where [the critical phase ~ k = t p - qo2] oscillates from those where it circulates. In each case ~k circulates in the region enclosed by the zone and librates in the region outside the zone. Trajectories follow a guiding trajectory [h°(I, ~, fi) = h] until they reach the zone of uncertainty where they exhibit some complicated motions until finally they reappear on some neighboring guiding trajectory (with a different action) and the process is repeated" [3c]. Fig. 3 (reprinted by permission) illustrates the beauty of Wisdom's explanation. The union of all averaged trajectories (from both regimes) which hit the uncertainty zone of energy h [h - UZ] forms a large ergodic region of motions with energy h. It is possible to transform the "principle" into a "theorem". Let T < ~ be the maximum of all times (in the z-scale) taken by averaged trajectories reaching RO)(ff, ~) = h from their "birth" to their "death". Fix a/J-neighborhood of that curve. Then there exists c = c(T,/J) and e0 = e0(T,/J) such that theorem 1 holds for e < e0, ,outside that neighborhood, which therefore defines the uncertainty zone. Precise estimates, however, seem hard to be obtained. Fortunately, as often happens in practice, the numerical results will most likely make sense beyond the range e0 guaranteed by such estimates. To make a distinction from another definition to be introduced in section 5, we may call this/J-neighborhood the "static" UZ.

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

98

It should be mentioned that for other resonances (such as 2/1 and 3/2, [18]), the instantaneous ¢p-system (25) (hereafter named the "frozen system" (FS)) contains more trigonometric terms and the transformation to action-angle variables (required for averaging) amounts to a tour-de-force of algebraic manipulations, to say the least. Nevertheless, we will show in the next section that most of the work can be circumvented if one is just interested in determining the extreme values of the eccentricity jumps. Before doing that we must state some trivial but useful remarks: i) Let q#l(x, y) the unstable equilibria of the FS, 1 _
yj,, l
R J l ( Y , y ) = R ( Y , fi, tpJ'(~,fi))=h.

(40)

ii) Action-angle variables are set around each stable equilibrium ~J2(x, y), 1 0, I = 0 at the center. Then the averaged Hamiltonian is of the form

h°(Y, .~, I) = R(~, ~, ¢pJ:(~, fi) + H(J')(Y, y, I), (41) H(J')(Y, y,O) -= O. The choice of signs used commonly in Celestial Mechanics make a < 0 so we must have H (j2) <_O. Hence the guiding trajectories of the j2-regime must be contained in the region

RJ2(~, ~) > h.

(42)

Unfortunately (42) is of very limited usefulness, in general. For instance, for the energy chosen in fig. 2, only a small region around (x2,0) is prohibited a priori. In contradistinction, the outer circulation trajectories (i.e., those corresponding to the FS curves external to all the separatrices) can be easily bounded. We illustrate with the 3/1 resonance. Write

h = H = F(x2 + y 2 ) + ffejx + h', Ix

(43)

h'= ½aA 2 + P c o s ~ . In the circulation regime h' < - P so

h < F(x2 + y2) + f f e j x - P = R(1)(x, y)

(44)

(we use unaveraged coordinates here). This is why the guiding trajectories of the circulation regimes are in the interior of the h-UZ central curve R = h. iii) From each point of R (1) = h emanates tangentially one averaged trajectory from each regime. The proof will be given in subsection 5.3 (proposition 4), but can be read immediately. iv) The curves h°(I, Y~,f,)= h (h fixed, I variable) cannot intersect, because h is a monotonous function of I for each (~, y) fixed. v) In practice, the boundary of the ergodic region can be obtained as the enveloping libration trajectory which ends at R (1) = h. The actual ergodic region is only slightly bigger.

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

99

4.2. Maximum eccentricity attained: 3/1 resonance Wisdom's diagrams (fig. 3, see also [3c, figs. 5-7]) are symmetric with respect to the x-axis. This of course follows from the fact that R ~ and P are even in y. The same will be true for other resonances, such as 2 / 1 and 3 / 2 [18], reflecting ultimately the symmetries in the geometrical configurations of the three bodies (see the figures in [6]). In view of the properties (v)-(vi) above, the enveloping libration emanates from the x-axis. Strictly speaking, one cannot rule out the possibility that the maximum eccentricity is attained at symmetric pairs with y ~: 0 but the diagrams indicate that it occurs at the x-axis, in the right side. Taking this for granted, we illustrate now how the extreme eccentricity can be roughly estimated. For concreteness, we choose as the total energy that of Sinclair's orbit Rl(xs,0) = -1.93 ej2 (see (19)). Substituting e j = 0 . 0 4 8 and multiplying by /~ we obtain H = -4.25 × 10 -6 which is very close to the energy value - 4 . 2 2 × 10 -6 considered by Wisdom.

Step 1. Using the harmonic oscillator approximation, determine the action of the libration trajectory trough (Xs, 0): I*

2 A(xs'0) =

1/2

(45a)

•

Step 2. Solve the algebraic equation in one variable R(2)(x,0)

-

2lCx 2 + D e j x + Ee211/21Cx2 + D e j x s + Ee~l 1/2 = R ( 1 ) ( X s , 0 )

(45b)

for a root x M > x s. Equivalently, solve S(u) = T(u), where

T(u) = - 1.93 + 2fllCu 2 + Du + El 1/2 S(u)= ( (F+C)u2+(ff+D)u+E'

uf~(ul'u2)' u ~ ( u l , u2),

(F-C)u2+(ff-D)u-E, u=

x/ej,

u 1---0.14,

(45c)

fl = ICu~ + DUs + E[ 1/2 = 1.26, u 2=2.93,

us=1.86

(see (14),(17)).

The table below summarizes elementary arithmetical calculations. We use simply the approximate relation e = ~/(l*la.) 1/a. There is a root between 0.3 and 0.4, as expected from Wisdom's numerical findings. Probably this is much too good, resulting from a lucky balance of two negligencies: omitting the transformation 2 ~ x and pushing the harmonic approximation for all the libration regime. u

e

S

T

5.2 6.9

0.3 0.4

5.5 14.7

6.0 10.3

We believe that closer agreement will be obtained using semi-analytical methods for action-angle

100

J. Koiller et aL / Relaxation-chaos phenomena in celestial mechanics I

reduction, such as [19], and performing the transformation (Y, y) ~ (x, y) via (O=0)

x=7,-eGy(x.,y), y =fi + eGx(Y, y ) .

See remark (ii) in the end of subsection 3.1, eqs. (28c, d). We may also estimate the libration amplitude at the maximum eccentricity value. From (32b) we get A~p = [ a / A ( x M , 0)11/4 2V~]*

=

2B1/2 [Cu 2 + Du M +

(45d)

El"

Substituting u M - 6 it results A~p - 1 rad - 60 deg. This agrees roughly with the deviations from the stable equilibrium q02(x, y) in [3c, figs. 1 and 4]. The bottomline is that the harmonic oscillator approximation may work reasonably well far beyond its range. 4.3. Maximum eccentricity attained: outline for the other reasonances We indicate what has to be done in the ( p + 1 ) / p resonances, and what is needed for precise calculations. We plan to attempt them on the subsequent paper dealing with numerical aspects. Following [18] the Hamiltonian has the following structure: H=½a~2+#R(x,y,q~),

a<0,

R = Co+ C1 cosq0 + C2 sinqo + C3cos2qo + C4 sin2ep, Co = B I ( X 2 + y 2 )

+ B2X, (46)

C 1 = B 3x + B4, C 2 = B3y , C, = B s ( x 2 _ y2) + B6 x ..[_B 7 ' C4 = 2Bsxy + Bty.

The variables have analogous meanings to the 3/1 resonance, but more trigonometric terms are present. Reduction of (25) to action-angle variables can still be done with elliptic functions, but the amount of algebraic manipulations needed suggests using MACSYMA or REDUCE, or better perhaps using semi-analytical methods [19]. Even the much simpler task of studying the constrained systems is not very easy. The condition R~ = 0 yields an equation of forth degree (1 - u2)(C2 - 4C3u) 2 = [ C l U - 2C4(1 - 2u2)] 2,

u = sin ¢p.

(47)

Solving it by radicals, one obtains closed form equations for the branches R (1) >_ R (2) > R (3) ~ R (4)

(48)

J. Koiller et aL~Relaxation-chaos phenomena in celestial mechanics I

101

a Fig. 4. Phase portraits of the frozen system for the (p + 1)/p resonance; y = 0. a) ]Cll < 14C31, C 3 > 0; b) ]Cll < 14C31, C 3 < 0.

(here R (1), R (2) are the stable, R (3), R (4) the unstable) Hopefully, however, only the values on the x-axis are needed to determine the maximum eccentricities. For y = 0 we have R = C o + C 1cos qo + C 3cos 2¢p, CO=B1Xz'4"Bz,

C 1=B3X'4-04,

(49) C 2 = B 5 X2"4-BrX'4-0 7 .

The local extrema are given by cp = 0, or,

(50a)

R = Co + CI + C3, R ~ = T- C 1 - 4C 3 (upper sign refers to cp = 0) and c o s ~0 = -

C1/4C3

R=C O

8C 3

Ifx/C31 <-4 )

(provided

C?

16C2-C 2 C 3,

R~o~ =

(50b)

4C 3

Depending on the values of the constants B and the x value different FS phase portraits are possible. See fig. 4. As in subsection 4.2 one may resort to the harmonic approximation for preliminary estimates:

h(J)(~,O, I ) = R(J)(£,O) -laP(J)(Y~,O)11/21,

j = 1,2,

(51)

P
Step 1. Given a value h for the total energy, find the x-axis limits of the U Z by solving R(3)(£,0) = h,

R(4)(x,0)

=

h (52)

( R (3) = R (4) in fig. 4a, R (3) > R (4) in fig. 4b).

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

102

Step 2. For each root ~tJ of (52) find the action value associated to trajectories leaving the UZ at YtJ. I*=

R(J)(xu'0) - k 11/2,

I,~P(J)(~o,0 )

J=1,2

(53)

(R (1) > R (2) in fig. 4a, R o) = R (2) in fig. 4b).

Step 3. For each I * as above, solve for Ix~l > [xuI the algebraic equation R(J)(Y, 0) - [ aP(Y)(Y, 0)[1/21 * = h.

(54)

Then the maximum eccentricity corresponds to the greatest of the [YMI-

Remark. The transformation of variables Y ~ x seems to be a nuisance dictated by the necessity of performing the near identity canonical transformation, in order to average out the critical angle. Perhaps a direct treatment in the original variables would be possible, using a two-time-scales method. A starting point, suggested by prof. J.B. Keller, is presented in appendix A.

5. The dynamical systems picture for relaxation chaos 5.1. The normally hyperbofic invariant manifold. Unstable periodic orbits Using the intermediate time scale o = et, we get an equivalent ("dual") form for (5), in which the cp-motions are the primary focus of attention:

d2q)

do 2

=

-aR~,

dx do

eRy,

dy = eR x. do

(55)

For e = 0 a normally hyperbolic structure is present (see [20] for background on dynamical systems) with invariant manifold M o : ( (x, y, ¢p m-~o(1)(x, y ) , A = 0 ) l P ( x , y ) * 0}.

(56)

The dynamics is trivial - that of the frozen systems - so M 0 has coinciding stable and unstable manifolds (of dimension three)

W= WS(Mo) = WU(Mo) = [,.J frozen pendula separatrices. (x,y)

The persistence theorem of hyperbolic structures [21] implies that for small enough e there is an unique invariant surface M~, e-close in the C 1 topology to Mo. It is still hyperbolic, but we will prove that the three dimensional local stable and unstable manifolds WloS,~U (M~) are no longer coincident but they intersect transversally. As it is well known from dynamical systems theory, this fact has far reaching consequences. Technically speaking, the persistence theorem demands a compactness assumption on M~ which some-

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

103

times is handled by modifying the flow outside the region of interest [22]. In our case the requirement is met via the energy integral. Indeed, we now show:

Theorem 4. M e consists of a family of unstable periodic solutions. Proof. As M~ is a graph over the (x, y) plane, we may write Me:

fp =

l~?(1)(X,

Differentiating and substituting into (la) it follows that manifold of constant energy H--/~h is easily obtained as

"lh,e:

(57)

y) + e~(X)(x, y) + e2qj(2)(x, y ) + . . . .

h = RO)(x, y) + ~__~[,~(1)t~ _Cpy(1)Rx)2+ 0(#2), 2a ~Vx --y

~(1)

~_.

0. The intersection ~h, e of M e with the

(58)

where R x and Ry are computed at (x, y, qo°)(x, y)). So 3'h,e is a closed curve which stays within a distance 0 ( # ) from the center R (1) = h of the UZ. Let h be a regular value of R 0) (i.e., we omit Sinclair's value). Then for/t small enough the velocity vector e ( - Ry, Rx) along 7h, e does not vanish, and is a periodic orbit, as asserted. []

Remarks. i) Theorem 4 implies that the "static" UZ width 8 should be taken of order higher than O(e) so as to fit Yh, e "comfortably" inside it. (Notice: we are using unaveraged coordinates here.) ii) The O(#2) term in (58) depends explicitly on qj2 and can be easily collected (appendix B). iii) Near Sinclair's point we have two periodic solutions "/h,e = ~1h,e U 3'~,e (1, r = left and right) such that W u (3,h,e)ch 1 - - W s ('/h,e) r 1 - - W u (~'h,e)r and W s (~/h,e)ch The dynamical consequences are of the same type than in the case of one component.

5.2. The Melnikov function In order to measure the gap between WlS(Me) and WlU(Me) we fix a (x, y) value and consider a line L~ in (cp, A) plane, perpendicular to the unperturbed pendulum separatrices (see fig. 5a) fp = (p(1)(X, y) T- 2 arctansech (3,o), A = T- 2V sech(~/o),

(59)

Ot

v 2= I e(x, y)I, S,U at o = a (the upper sign refers to the upper separatrix A > 0). L~ will meet Wlo ~ (Me) at the points which we call p~.U(a, x, y) respectively. The distance between these points is measured in terms of the level values of

h" = ½aA 2 + Pcos +,

~ = ~p- qg(2)(x, y ) ,

(60)

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

104

namely

A,( a, x, y) = (grad(A, ~)h'lp ~ _ pS ) = eM( a, x, y) + 0(e2).

(61)

The function M(a, x, y) is called the (first) Melnikov function. Its importance stems from the fact that if grad(,, x, y ) M ~ 0 at a point where M = 0 then W~(M~) and WU(M~) must have to intersect transversally on a point near (x, y, ep(a), A(a)). In appendix C we compute M using C. Robinson's formula [22]. We show that M is independent of a, and in fact, for both separatrices, 4

M(x, y ) = Irl'a"'a/2 (Rsec, P ) ,

(62)

where ( f , g) =fxgy-fygx is the usual Poisson bracket. Since Rsec and P are even in y, M(x, y) is odd with respect to y. It follows immediately that the stable and unstable manifolds intersect transversally along the x-axis. More explicitly,

M(x, y) -~ 2yN(x, y)ej/lall/2p 3/2, (63)

N(x, y ) = 4 C ( x 2 + y 2 ) ( - F D + ffC) + 4ejxC(Dff-4EF) + 2e ] [ff(D 2 - 2EC) - 2FDE]. It follows that M also vanishes on a circle. Substituting the numerical values, we get

M( x, y) = -2.9859ejCyE(x, y )/lall/2P 3/2, E(x, y) = (x + 0.3083eJ) 2 + y 2 _ 1.2742ej2.

(64)

Remarks. i) The sign of M has an important dynamical interpretation (fig. 5b). When M > 0 trajectories may enter the libration regime, while when M < 0 may leave. ii) We must also explain why M should be independent of a. Recall that whenever M = 0 there is a nearby homoclinic (doubly asymptotic) orbit to 7h, ~. Since thevelocities dx/do, d y / d o are d~(e), the first Melnikov function seems to be unable to detect the initial changes in x, y as a varies along the separatrix. An algebraic explanation is given in appendix C. We presume that the second Melnikov function will depend on a, but its calculation is not necessary for our purposes here. iii) For ej = 0, M vanishes identically and one expects the system to be integrable. This is indeed the case as we will see in subsection 5.4. Notice that the condition M :~ 0 is equivalent to Wisdom's criterion ( R s¢c, P } #= 0 for chaos.

5.3. Dynamical systems interpretation for the origin of chaotic behavior According to dynamical systems theory [20], transversal intersection of WS(Th,~) and WU('th.~) along homoclinic trajectories implies the existence of random motions and nonintegrability of the system [23]. In frequent applications of the Melnikov method (see e.g. [24]) one considers small periodic perturbations of integrable systems, H = H°(p, q) + eHl(p, q, t), and one studies Poincarr's one-period mapping. As it is well known, a thin chaotic region appears around the separatrices in the (p, q) plane (say, the phase space of a pendulum (% A)).

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

L 8

M>O

b

M
Fig. 5. a) Melnikov function ingredients, b) Change of regimes (schematic, (x, y) drift not represented).

0.2

0.1

y 0.0

I

I

l

-0.1

-0.2 -0.1

I 0.0

I

I

0.1 x

0.2

0.3

Fig. 6. Projection in (x, y) plane of a stretch from an ergodic trajectory (reprinted from [3c, fig. 8] with permission).

105

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

106

In contradistinction, in the present context it is more important to project the four dimensional trajectories with constant energy in the (x, y) plane of the "slowly varying parameters". As Wisdom has emphasized, the most important feature of the relaxation chaos mechanism is that these projected ergodic regions in the (x, y) plane are quite extense. From the geometrical viewpoint, what happens is that the global stable and unstable manifolds (which are Lagrangian submanifolds) form a complicated web of intersections projecting densely over a (x, y) region. The chaotic trajectories are those trapped in this "labyrinth" inside the three-dimensional space of constant energy (fig. 6). The "classical uncertainty principle" asserts that the dynamics of any solution of (1) inside the ergodic region can be divided, for the sake of numerical investigation, into clear libration and clear circulation periods. The transition motions form the UZ, which thus plays a role analogous to the "boundary layer" in singular perturbation problems. The uncertainty Zone width 8 was itself left uncertain in [3c]. Recall from subsection 4.1 that the constant c in theorem 1 grows as ~ - , 0 so the appropriate choice of ~ seemed, at that point, just a question of experimentation. On the other hand, the dynamical systems approach suggest that another, a "dynamic definition" of the UZ could be provided. We can at least give a lower bound for 8. Our ansatz is based in the observation that, since the full ergodic region reflects properties of the global stable and unstable manifolds Wg~'o~(Yh,e) , then by the same token, the UZ should reflect properties of Wl~U(yh,~), and so should be taken at least of order 0(e). In fact, using the dynamical meaning of the Melnikov function, we propose the following

Definition. The (primary) uncertainty zone of energy h is the annular domain Ih - R(1)(x,

y)[< elM(x,

y)[.

(65a)

So the width of the primary UZ at a point (x, y) in the center curve R (1~= h is (65b)

d = l M ( x, Y)el/llgradR(1)( x, Y)II.

Note that inside the primary UZ (65a) the averaged trajectories cannot be safely used to approximate the motion. Indeed, we denote as in (60) h' = ½aA: + Pcos ~p. Then h - R (1) = h'

+ Rse c -

R O) --- h'+P.

(65c)

Thus IR(1) - hi < Me is equivalent to [h' - ( - P ) [ <

Me.

(65d)

Suppose that p~,S are in opposite sides of the unperturbed separatrix. These inequalities imply that (¢p, A) is in the segment between p~ and p~ where trajectories are captured for the other regime (see fig. 5a). A comment is in order here. The trajectory through p~(a, x, y) tends as t ~ + oo to the periodic orbit 7~(h(pS(a, x, y)) whereas the trajectory through p~(a, x, y) tends as t ~ - oo to 7~(h(p~(a, x, y)). These

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

107

energies are in general different. We suspect that this line of thought can lead to the proper "dynamical definition" of the h-UZ, but so far we have not been able to frame it properly. We cannot refrain ourselves from making a preliminary guess, though: the dynamical h-UZ boundaries consist of all (x, y) for which there is a parameter value a such that (x, y, pU'S(a, x, y)) has energy h. We observe that in the "regular" case of homoclinic behavior, stable and unstable manifolds, together with the Melnikov method, have been successfully used to determine the boundaries of the thin ergodic region [25]. We suspect that the present situation is similar, although more involved. We also wish to call attention that, outside the extense ergodic regions in the (x, y) plane produced by the relaxation chaos mechanism, the numerically obtained averaged trajectories look as regular closed curves. For the unaveraged motion this should unfold into invariant toil, in which the frequency of ~-motion is 0 ( l / e ) times that of (x, y). Is then the case that KAM theorem (i.e., existence of invariant tori) may hold in some situations where the nondegeneracy hypothesis is not satisfied, and worse, in which the frequencies ratio depends on the perturbing parameter? In connection to some geometrical reasonings we will present in section 6, we now prove a simple result which was mentioned earlier.

Proposition 4. The averaged motions h = h°(I, ~, ~) approach tangentially the center curves R O) = h. Proof. The ODEs for averaged motion can be written as dx

d---~= - 2Fy - Py(cos ~k), (66)

dy d--~ = 2Fx + Fej + P~(cos ~k),

where (cos~k) is the average value along the instantaneous pendulum trajectory (25) with action I. Following [3c], 2E(kL) K(kL ) (cos ~k) =

1,

for librations,

2E(kc) 2 k 2 K ( k c ) + 1 - k--~c, for circulations ,

(67a)

where E, K are respectively the complete elliptic integrals of the first and second kinds (see [26]). The elliptic modulus satisfies R (1) -

k2L=

h + 2P

2P

'

k2

R (1) -

2P h + 2P'

(67b)

and moreover ip ig = 8_ -a'71"

1/2

ic=41P ~-

x/zg(kc) k---~'

[E(kL)-(1-k2)K(kr)]

'

~.laPi1/2 ~0g= 2 K ( k g ) (67c)

~rl(RO)-h+2P)a/211/2 ~°c= K(k¢)

108

J. Koiller et a l . / Relaxation-chaos phenomena in celestial mechanics I

In reality R (1) = h, the the average contribution ODEs.

these formulas are not needed immediately. Just notice that as (ff(~-),)7(~-)) approaches elliptic modulus (in both regimes) approaches 1. Since E converges and K diverges as k ~ 1, value (cos ~ k ) ~ - 1. Intuitively, for the "average" of cos ~k at the separatrix only the at the unstable point is necessary. Now, if we set (cos + ) = - 1 in (66) we get branch 1 []

Remarks.

i) It is easy to see that trajectories of R (1) run clockwise (see section 2). Hence the averaged trajectories in fig. 3 must run from y < 0 to y > 0 in the libration case and from y > 0 to y < 0 for circulations. This is consistent with the sign of Melnikov's function (fig. 5) in case the UZ is crossed, ii) Figs. 6 and 12 of [3c] seem to be conflicting near R o) = h, but we believe this is due to the nondifferentiability Of the mapping (~,)7) ~ (x, y) at R (1) = h: the "corrected" fig. 12 is a distortion of fig. 6. 5.4. Eccentricity of primary taken as the small parameter In some other situations of resonance (e.g., the motion of a planet in the gravitational field of a double star) it may be more natural to consider/~ a fixed quantity and ej as the small parameter. Here we will touch this case very briefly, only as a means of comparison with the relaxation chaos mechanism. The Hamiltonian is now given by H= H ° + ejn I + 0(e2),

(68)

where H 0 = ½a~2 + i.tF(x 2 + y2) _ t x C [ ( x 2 _ y 2 ) c o s ~

+ 2xy sin ~0],

(69)

H 1 = # fix - I~D x cos q~ -/~ Dy sin ¢p.

The Melnikov function in subsection 5.2 suggests that H ° is integrable. Indeed, using polar coordinates (x, y ) = 2 ~ ( c o s w, sin w) we have H ° = ½acb 2 + # r z p - # C 2 p cos (2w - ~)

(70)

and the canonical change of coordinates

0a=2w-(p, 02 =

2W,

11=-~, (71)

I2 = dp + p / 2 ,

yields

(

H ° = ½~ I1 + T

- 4 1 ~ C ( I 2 + I 1 ) c ° s O 1 + 4F/~I2

so that 12 is a constant motion.

8F2/~2~,

(72)

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

109

For the perturbed system one expects the occurrence of the regular type of homoclinic behavior. In fact, the simplest among the equivalent Melnikov functions for this problem can be written as (see [27])

U(a) =

(om/

I73)

O2)d,,

where the integrand is evaluated along a separatrix ( 0 1 ( t - a), i l ( t - a), 0 2 ( t - a)) of the unperturbed system with a fixed value of 12. Let us just indicate that actual computation of (73) by residues should present no difficulties, since the separatrices can be written in closed form in terms of elementary functions. To simplify notation, introduce P = 11 + 4F/~

4F/~

Ot

Y = I2 -- - Ot

q

----- 0 1 '

a/3 = 4/~C,

(74)

so that H ° = ½ap 2 - aft( p + y) cos q.

(75)

Thus 1 ~-0 = P - / 3 cos q.

(76)

Solving (75) for p with H ° --- ah and inserting into (76) we get +_a/3( t -- to)

(

dq . / ~ ( c o s q - y//3)2 + 2h - y2//32

(77)

The substitution u = cos (q/2) yields du . f ~/(1-.2)[(2u2-1-v//3)2 + 2h-yz/flz]"

_+ - ~ ( t - t o ) =

(78)

Now, the equilibrium points of (75) are (i) ( p = - V , c o s q = - 7 / / 3 ) p r o v i d e d Ivl < 1/31; (ii) ( p = fl, q = 0);

(79)

(iii) ( p = -/3, q = ~r). One verifies easily that if 171 < 1/31 then the solutions (i) are unstable, while (ii, iii) are stable, and if IYI > Ifll then one of (ii, iii) is unstable, the other stable. For instance, the unstable orbits (i) correspond to h = ½"r2 so (78) becomes +_afl(t-to)=

f

u2

du l + y / f l ( 1 - u 2~1/2"j 2

which is an elementary integral.

(80)

110

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

a

!

s!

(x,y)

b

(cO, A) Fig. 7. An averaged "string" of motions looses up upon hitting the UZ.

6. Discussion: how to simulate motions inside the UZ?

The ideas presented in this section are very tentative. Our aim is solely to highlight the features which bring about the global chaotic behavior. Numerical experimentation and/or revisions will be reported in the subsequent paper. The basic idea we have in mind is eventually being able to simulate UZ motions by "random composition of mappings", a mathematical topic which has aroused considerable interest in the last few years (see e.g., [28]). A different approach might be possible through the method outlined in appendix D. To fix ideas, suppose that a trajectory enters the primary h-UZ at the point A (fig. 7a). Of course, no information is known about its phase O (which was averaged out). Unfortunately, the exit point B should be very sensitive with respect to this phase - ultimately, this must be the source of randomness: the would be invariant tori which "scatter" upon hitting the UZ. (fig. 7b). Notice that even if the phase at B were "measured" (say, by astronomical observations) it would be lost during path BA' (fig. 7a). In fact, the variation of 0 along BA' is

A0 = ~-

a,W d r + JBA' aI d r + 0 ( 0 .

(81)

Informally speaking, not only 0 will wrap many times around S 1 but also the fluctuations due to the perturbative terms will become important, so it is impossible to keep track of the phase. For an anecdote

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

111

reminiscent of this situation, we quote Feynman [29, "lucky numbers"]: "one day I was feeling my oats [and] I announced, 'I can work out in 60 seconds the answer to any problem that anybody can state in 10 seconds, to 10 percent' [At that moment], Paul Olum walked by in the hall ( . . . ) . Without hardly stopping, he says, 'the tangent of 10 to the 100th'. I was sunk: you have to divide by pi to 100 decimal places. It was hopeless". Let K be the set of all allowable libration or circulation actions (i.e., of those averaged trajectories reaching the UZ, K = K L k.J Kc). We will now define a mapping F0o: K ~ K, F0o(IA) = action at exit point B,

(82)

where 0o is the phase at entrance point A (a random variable on $1). The exit point B(A, 0o) is obtained in two steps: (i) regressing to the standard phase, which corresponds to the "lowest" pendulum position; (ii) constructing a map U for successive standard phase points. Then B will be the first iterate emerging out the UZ. [One of the reasons why B is very sensitive to 00 is that the number of iterations jumps at certain values of 0o.]

Step 1.

Regressing to the standard phase. From (26) it follows that the standard phase is 0 = 0 for the circulation and 0 = 0 or ~r for the libration regime. For definiteness suppose that we are arriving at A from librations. During a short time there is no harm done in using the averaged equations d£ o d'-'~ = - h Y '

dfi _ o ~-h~,

dO l h ° ( ~, fi, IA) d--~=

(83)

backwards in time, so we want the first TO< 0 such that

1{~ow Tj °

[ -0o, dT=[_0o+~r '

i f 0 < 0 o<~r, if~r<0o<2~r.

(84)

Let (x o, Y0) = (x(T0), .V(To)). Then we have

Proposition

5. a) 0 > To > -~re/o~(A), ~o(A) = const le log eI.

~rlap[1/2/2K'([Me/2P[1/2);b)

IR(1)(Xo, Yo) - R(1)(A)[ <

Proof. Since

~o",~0 as we approach R 0) = h, we have ~o(T) > ~0(A). Then a) follows immediately. For b) we give a more precise estimate

IR(1)(x°'Y°)--R(1)(A)]<

[

Jo °

gradR(X)"

d r ' dT

dT < [To[ max -

grad R (1)- {

o<~<1~oll

Now, (66) implies

dT'

"(-Py, Px)(1 + cos ) [ 21(Rsoc, P) I

d~

~ dT

d Y\ ] ' dT

j

"

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

112

Fig. 8. Coordinates for the UZ motion mapping.

and hence max I{ R~c,

IR(1)(Xo, Yo) - R(1)(A) [ <

P ) l'4eg'(IMel2ema~l) [~pmm[1/2

'

(85)

where Pm~, = maxR(~=hP, Pmi, rl~Rtl)~hP and K ' is the associated elliptic integral of the second kind. Recall the classic formula K ' ( k ' ) - log ( 4 / k ' ) [26]. =

Proposition 5 suggests the following (still tentative) Definition. The secondary uncertainty zone (IIUZ) is the annular domain Ih - R ' ( x , y ) l < mhe + cheK'[( dhe)l/2],

m h = max IMI, R(1)=h

dh=

c h = 4 max [{Rsec, P } I/laem~nl1/2,

(86)

RO)~ h

m J 2 P m ~ ,.

So the secondary UZ width is O(elog e) as e ~ 0. Compare [12b, §5.5]. Remark. Proposition 5 shows that the regression path (~(A, 0), .~(A, 0)) has time length A~" and "transversal variation" Av = R
~1 = % + A¢.

(87)

J. Koiller et al. / Relaxation- chaos phenomena in celestial mechanics I

113

Starting from the unaveraged equations dx/dr

= - 2 F y - Pycos ~p- Pcp~2) sin +,

d y/d,

= 2 F x + f f e j + Pxcos ~ + p~(2) sin q~,

(88)

consider one non-periodic oscillation of ~o near the separatrix (i.e., one circulation or one libration swing). Proposition 5 shows that the variation of (x, y) is only of order O(elog e) so the contributions of the sine terms tend to cancel by symmetry. We drop these terms so that d v / d ~ ' = grad g 0).

dx d y ) dr' ~ = - (Rse¢' P } ( I + cosqJ).

(89)

Our idea its to divide the pendulum oscillation into two parts: the "upper", where cos ~ < 0 and is slow at the top (so we may approximate roughly cos ~k by - 1 and hence Av = 0); and the "lower", where it is quickest at the standard phase point (x 0' Yo) and where Av is to be produced. According to the Machiavellian principle that "all evil should be done at once" we concentrate the variation of (89) at (x 0, Yo) and we approximate the lower qJ oscillation by a piece of the FP separatrix. An elementary calculation gives fc

os4,>O

( 1 + cos qJ) do =

(on the separatrix).

7

(90)

It follows that Av= -e2v~

(Rse¢'P}

-

(91)

eM(r°)

where M(%) denotes the Melnikov function evaluated along the central curve. For A,r we take the average time of the frozen pendula oscillations in the interval [Vo, vo + Ao]. Hence, from h - R (1) --- h' - ( - P ) we obtain 2~ d o ] - ]12c°s ~d = y2(1 - f l ) '

(92)

fl -- IvIP'

SO

v~- f~os-kl-I,BI) T(v) = ~

"0

d4'

(cos

q~+ (1

= 2K,(~). -1/31)

(93)

Thus it follows that (this formula is valid for both regimes)

laP[ Av

)

do.

(94)

Remarks.

i) As expected, A~- grows near the central curve o = 0. That is why the exit point becomes so

114

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

unpredictable. Whenever M(r0) = 0 then we may set

Ar = 2eg'( (IoV2p )l/2)/lae I.

(95)

If v = 0 also, then Az = oo and U becomes undefined. Perhaps this is a reflection of the fact that in numerical experiments (see e.g. Figs. 1-4 in [3c]) the eccentricity apparently falls to an apparently random value after reaching its peak. We suggest using a random number generator for Ar whenever U is undefined. ii) Numerical evaluation of (94) should present no difficulties. If oo, o I are in the same side, the trapezoid rule and a table of elliptic functions are sufficient. If v0, o 1 belong to opposite sides it requires a little more care. Use Ar = ~_~_~ o 2 e {[ j flool^ . . . .aVtoflOllK'dv) .

(96)

and a trunction of the well known expression (so that a closed form can be obtained for the integral) K'(k')=

log4+

(12 )~2/[log~-7 4

7 .2 ) k ' 2 + / 1[2-~) " 3 x2/[l°gk4"7

2 -4-~.4 ) kt4 1.2

Jr'''

(97)

iii) For e infinitesimally small and opposite vo, v 1 we may retain only the first term in (97) so e Ivollloglv01l +tvlllloglvlll At- -- lael Ivol + Ivll

(98)

iv) Since each iterate of the mapping changes v by O(e) and the IIUZ width is O(e log e) we expect that the number of iterates needed to cross it to be proportional to log e. Roughly,

(m) = ~

'l'=h I M ( ' ) I

loge.

(99a)

We estimate the average iteration time by

(AT)=

arealIUZ

,uz~¢(¢,o)d~'dv--- ~-~ ~

-~- dlogd.

(99b)

Hence the average time of permanence of a trajectory inside the IIUZ is 2v~-c h

,o,

1

1

dr

2

,00)

Heuristically, for e = 10 - u the number of iterations grows proportionally to N but the UZ is crossed in time of order N 2 X 10 -N.

J. Koilleret al./ Relaxation-chaosphenomenain celestialmechanicsI

115

Acknowledgements We thank SBMAC (Brazilian Society for Applied and Computational Mathematics) for making possible communication among the authors, living in different cities, during the first stage of the work in 1985. Encouragement and/or sympathetic raised eyebrows were given at several occasions by Profs. Silvio F. Mello, Sueli Guillens, Roberto V. Martins, J.L. Sagnier, W. Sessin, Hildeberto Cabral, and Sonia P. Carvalho. One of us (JK) wants to thank Prof. J.B. Keller who read a working draft of the paper and suggested a two-time scales approach for averaging; Prof. P.J. Holmes who gave us the idea of exploiting a transverse hyperbolic structure; Prof. J.E. Marsden who informed us about the status of the exponentially small horseshoes problem; Prof. J. Wisdom who explained us his approach in a sunny Sunday afternoon; and Prof. Vincent Moncrief who sponsored our visit to Yale University Mathematics Department. Last, but not least, to the anonymous referee who suggested using Whittaker reduction in further developments.

Appendix A Outline of a two-time scales method The following is a suggestion of Prof. J.B. Keller, based upon [30]. Recall the equations of motion in the "r-scale

xr= - R y ,

y,=Rx,

Define e = IX1/2, O =

e-l'r,

x=X(,r,o,e),

g%r = - a R ~ . and let

y = r('r,o,e),

=

o,

(A.a)

Then

X,+e-IX,,=-Ry,

Y,+e-IY,,=Rx,

e2(~,~+2e-'~,o+e-ztboo)=-aR~.

(A.2)

Now let X = X°('r, o) + eX°)(T, o) + O(ez), etc.

(A.3)

Substituting in (A.2) it follows that

X°o=O,

Y°=O,

dP~°=-aR,(X°,~/°,~P°).

(A.4)

Hence X ° = X°(~-), y0 .= y0('r) and

o) =

_

+0('r, o)).

(A.5)

The general solution of (A5) has two integration "constants", a phase ~('r) and an energy E('r). Thus , o = V [ o + ~('r), E('r); X°('r), yO('r)].

(A.6)

116

J. KoiUer et al./ Relaxation-chaos phenomena in celestial mechanics I

We assume that V is periodic in o with period P(z). Now the term of next order is e and (A.2) yields:

X(o'= - X ° - Ry[X°,Y°, O°],

(A.7)

Y,('= - Y° + Rx[X°,Y°, ~° ]

(I)(1o)+ aR~[X o, yO, ~o ] (/)(1)= _2(/)o0 _ aRx~(X o, yO, ,~o)X(1)_ aRy~(X o, yO, ~0) y(1).

(i.8)

We require that X (1), y(a), ~(1) be periodic in o with the same period P(~-) as ~0. Then integration of (A.7) over one period P ( z ) in o yields

X° + (1/P('r)) fo

P(r)

(A.9)

Ry[X°('r),Y°(T),cb°('r,o)] do=O,

yO _ (1/p(,))foP(*)Rx[XO(z), yO(~.), 0o(~., a)] do = O.

(A.10)

Remark. Notice that equations (A.6) and (A.9, A.10) are still coupled through the energy function E(,r). In hindsight, the uncoupling must come from the "adiabatic principle", but it should be interesting to carry out the next terms in the expansion to see if it comes also from the formalism. We believe, however, that the Hamiltonian character of the equations must be explicitly taken into account.

Appendix B

The center manifold We substitute (57) in the equation %° = - a R ~ . The right-hand side is

-ctR~ = - a [ - P ( e ~ (1)+ e2~ (2)+ " " ) + 0(e3)].

(B.1)

The left-hand side is obtained by differentiating (57) twice: q)o = 8[--fp(xl)Ry --I-ep(yX)Rx+ e(-q~)Ry + q/yl)Rx) + e2 ( - ~p(~2)R,+ ff(~2)Rx) + . - .

],

(B.2)

%o= e2S + . . . , S = r°(X)R2rxx--y-It-cp~yRx(1)2 _ 2eP(x~RxRy _ Cp(xl)[_RyxRy+RyyRx+Ry~(_Cp(xl)Ry+q)(1)Rx)] + ~ ( y l ) [ - R x x R y nt- R x y R x q- R x ~ ( - ~ ( 1 ) R y q- ~(1)Rx) ] .

(B.3)

Equating powers in e we obtain ~b(1) = 0,

q~(2)___S/aP,

(B.4)

where qo is replaced by ~(1)(x, y). Using now the energy integral, we obtain (58). To collect the d)(# 2) term, observe first that ~9 = (--q?(xX)Ry - ~(yl'Rx) ~'t-~2[(-~(l,Ry~--t- qo(l'Rxqo)@ (2,- @(2'Ry + ip(2)Rx] + -..,

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

117

where again ep = cp(1)(x, y) in all derivatives of R. Hence the/~2 term has a contribution

IOl (\ - - Wx ( ' ) R y .=~=(~(yl)Rx ) [ ( -~(xl)RyTo -[- (~(yl)Rx~ ) +(2) - ~)(x2)Ry -.~ I/)(y2)Rx]

(B.5)

arising from the "kinetic energy" and a contribution P(~b(2))2/2. arising from the "potential energy" - P cos (/~ q~(2)+ . . . ). Notice that (58) is a series in/~, not in e.

Appendix C

Melnikoo ' s integral We recall C. Robinson formula [22]. Consider a system of ODEs x ' = f ° ( x , w) + efl(x, w) + 0(e2),

w'= er(x, w) + O(eZ),

(c.1)

where x ~ R 2, w ~ R N, f o = ( - H q, l i p ) is Hamiltonian, H = H ( p , q, w ). Assume that the unperturbed system possesses a family of separatrices if(t, w), lim Y(t, w) = x~(w). t-..* + oc~ Hence there is a normally hyperbolic structure for e = 0 whose N dimensional invariant manifold M 0 = {(x~(w), w)lW~ R N } has coinciding stable and unstable manifolds W = WU(Mo) = W~(Mo) of codimension one in R 2 x R u. Let M~ the perturbed invariant manifold and WS'U(M~) the stable/unstable manifolds. Fix a w-value and consider a line segment L~ in x-plane centered at Y(a, w) perpendicular to f°(~(a, w), w). L a will meet W su ' ( m ~) at points x~'~u(a, w). One measures the separation between these points by A ( a , w) = (grad xHIx u - x~) =eM(a, w) + O(e2),

(c.2)

where M is given by

M ( a , w ) = f f % ( g r a d x H I f l + %f--~'v(t))dt.

(c.3)

Here grad~ H, f l and the 2 x N matrix Of°/Ow are computed along the separatrix trajectory ~(t, w) and

v(t) = fatr(Y~(S, w), w) ds.

(C.4)

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

118

Notice the useful fact that

M(a,w)= M(O,w)+ [f~_%(grad~Hlafo/aw)dt]. foar(Y(s,w),w)ds.

(c.5)

Melnikov's criterion asserts that if M has simple zeros (i.e., grad M 4:0 whenever M = O) then W~s and W~u intersect transversally at nearby points. In many cases f~oo(gradx H[Ofo/ao~)dt= 0 so that M will depend only on w, not on a. Some of these situations will be given at the end. Calculation of Melnikov's integral for the 3/1 resonance: x=(A,~0),

o~=(x,y),

H=½aA2+Pcos(q~-c;(2)(x,y))+Rsec,

(aA,-Psin~), ~p=cp-(p(Z)(x, y ) , f°=(Psin+,aA), fl-O, r=(-Ry Rx)= ( -2Fy- PyC°S+- Psin~kq)(2) ) ' 2Fx + ffej + Pxcos + + P sin +~o~(2) ' Of° (Pxsin~b-Pcostpcp(~2) Pysintk-Pcos~pCp(y2)) gradxH=

Ow=

0

cos cp(2) =

-A/P,

0 sin q)(2) = _

cos 1/, = 2seeh 2 yt - 1,

(C.6)

'

B/P,

sin ~k= T-Zsech(yt)tgh(yt)

(upper signs refer to the "upper separatrix" of the unperturbed pendulum) Robinson's formula yields ½M(a, w) = f_~% • sech u ([ Px( T-2sech utgh u) - P(2 sech2 u - 1)cp(~~)] f ( u )

+ [Py(T-2sech u tgh u) -

P(2 sech2 u - 1)~py(2)] g( u)} du,

(C.7)

where f ( u ) = Lu/v[ -

2Fy -

Py (2 sech2 7s - 1) - P(-T- 2sech ys tghys)Cpy(2)] ds,

g(u) = Lu/V[2Fx+ ffej + Px(2 sech2 y s -

(c.8) 1) + P ( • 2 sech 3,s tghys)cp~(2)] d s .

All integrals involved are elementary. The manipulations are alleviated observing that we are in one of the cases where M independs of a. Indeed, disregarding the odd integrals in the factor inside brackets in (C.5), we end up with the computation of f_~ sech u(2 sech2 u - 1) du = 0. 0C

(C.9)

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

119

Calculation of f ( u ) and g(u). Let v = ys. Then

f(u) = ¥1 rrj.l_ 2 F y g(u) =

fU[2r

Py(2 sech2 v - 1)_+ 2P~(yz) sech v tgh o] dr,

+ Fej + e (2 sech

-

-T-2 Pep~2)sech v tgh v ] dr.

Hence

y f ( u ) = ( - 2 Fy + Py ) u - 2 Py tgh u -T-2 P ~ 2) sech u, "tg( u ) = (2Fx + f f e j - Px) u q- 2P~ tgh u + 2Pqo(xZ)sechu

(C.10)

(the constant terms are irrelevant due to (C.9). Substituting into (C.7) are retaining only the even integrands, one obtains, after some nice cancellations, the desired formula (62). Some situations in which M independs on a: i) H = p 2 / 2 + f ( q , w), where f is even in q. Reason: taking the origin of time in the separatrix at the symmetry point, then the integrand in the expression inside the brackets in (C.5) is odd. ii) H = p 2 / 2 + A ( w ) c o s ( q - q(w)) (this is our case here; more generally, this can be thought as a particle moving on a wave with slowly varying amplitude and phase velocity). Reason: the question boils down into showing that f~p p c o s ( ~ k ) d t = 0 , where ~k=q - q(w). This either follows from the "cosh miracle" (C.9), or from the observation that for f = c o s ( q - q(w)) we are just calculating fpfqqdt = f~ dt = [P]-~oo = 0. (This explains the "miracle"). iii) The homoclinic loops in (p, q) plane are of "figure 8" type, and symmetric with respect to q-axis. Reason: origin of time can be made at the symmetry point. Then p(t) is odd, while q(t) even. The integrand Pfqw is therefore odd. We ignore whether there are cases in which the first Melnikov function depends on a.

Appendix D Whittaker reduction The referee has kindly pointed out that (1, 2) may be rewritten

H=~a~b2-l~Ee2cos~+tt[F-Ccos(2o~-~p)]2p+#ej[ffcos~-Dcos(~-~p)]V/~

(D.1)

(where we recall p - ½(Pla.)l/2eZ). The eccentricity being the variable of interest, it is natural to solve for p = p(H; o~, ~, ~), which satisfies a quadratic equation. Thus H - a ~ 2 / 2 + #Ee 2 cos ~o e2 [ffcos o~ - D c o s (to - ep)] 2 + 2 0 = b t [ F - Ccos(2o~-cp)] 2 I F - Ccos (2o~ - q~)]2 +

pEe 2 cos ¢p + Ccos (2~ - ~0)1

H - a~2/2

p[F-

+

[Fcos

4-~-

Dcos

Cc~~--)~

- , ) ] 2 1/2

ej f f cDcos(to~).ff_L___C_c~ss( oso~~-~--

.

(D.2)

J. Koiller et al. / Relaxation-chaos phenomena in celestial mechanics I

120

The referee suggested examining Whittaker's method [31], which allows the reduction of any n degrees of freedom autonomous Hamiltonian system to a nonautonomous n - 1 degrees of freedom system. Here we get dep d~o

ap 3~'

d~ d~0

0p ~¢p"

(D.3)

If one were able to solve (D.3), either analytically or numerically, for initial conditions ¢Po,~o, Wo, the resulting functions = cP(CPo, ~o, °~0; w),

~/i = q~(CPo, q~o, ~°o; ¢)

(D.4)

could be inserted into (D.2), so that e2 = e2(60)

2 p(H; ~0, ~0(¢po, q~o, ¢°0; w), qli(¢po, ~o, °~o; w)). (l~la.) 1/2

(D.5)

Let us therefore examine (D.3) in more detail. The situation is best understood if we rescale variable to A = ~ / e , and the energy equation to h= H

aA 2

EeZcoscp+ [ r - C c o s ( 2 o ~ - c p ) ] 2 p + e j [ f f c o s o ~ - D c o s ( o ~ - ~ p ) ] ~ / ~ .

(D.6)

Whittaker reduction modifies to

dA do~

1 Op -~ - ~ ( ~ , A, w; h) =

dcp do~

1 3p e

- ~-(¢p,

A

'

o~; h ) -

1 Oh/Oep e Oh/Op ' 1 ~h/~q~ e ~h/~p

-

.

(D.7)

(D.8)

Hence, if we rescale o~ to ~ = ~0/e (a "slow angle" (!?)) we get dA 3p (ep, d~ = ~ A,e~;h),

dcp 3p ~=,-~-(¢p,A,e~;h).

(D.9)

These equations have the nice feature of being cast in the usual form for the adiabatic principle. As long as one is away from the separatrices of the frozen systems P(% A, ~0; h), good approximate solutions can be obtained by transforming to action angle variables (J, a) and postulating the adiabatic invariance of J. As before, the method breaks down if, as the system evolves, the area inside the separatrix approaches 2~rJ. Equations (D.9) have the further advantage that no coordinate transformation was performed (compare subsection 3.1). Unfortunately, however, obtaining J(cp, A) seems to be a difficult exercise on special functions. Notice that since

Ee 2 cos cp eZ[ffcosw - Dcos(o~- ~)] 2 F- Ccos(2o~-~) + 2 [ F - Ccos (2o~ - cp)] 2

h - aA2/2 +

2p= +_

h - aA2/2 + Ee~ cos ep e ~ [- f~f-cCo -scwo-- ~D- -c-o~si(~o ~ - ¢p)] 2 1/2 e, f f c o s w - D c o s ( w ~ ) -P~-Cc~s-Ci-dF - Ccos (2w - cp) +

, (D.10)

J. Koiller et aL/ Relaxation-chaos phenomena in celestial mechanics I

121

t h e n A 2 satisfies a q u a d r a t i c e q u a t i o n whose coefficients c o n t a i n r a t i o n a l functions of cos cp of degree four. O n e is c h a l l e n g e d to calculate a n a l y t i c a l l y the action

W e also w a n t to m e n t i o n that there is a m i n o r technical difficulty, n a m e l y t h a t W h i t t a k e r r e d u c t i o n c a n n o t b e a p p l i e d at the turning p o i n t s where d w / d t = 0, i.e., when Oh/Op = 0. W i s d o m d i a g r a m s show t h a t this will o c c u r for a large class of initial conditions. I n such cases the range of w = e~ will b e a s y m m e t r i c interval. This p r o b l e m can b e c i r c u m v e n t e d b y m o v i n g a l t e r n a t e l y f r o m one to the other r o o t (D.10) a n d i n t e g r a t i n g f o r w a r d s or b a c k w a r d s . O n e s h o u l d n o t b e m i s l e a d into t h i n k i n g that this is a cause for the e c c e n t r i c i t y increases; rather, this is j u s t a h a r m l e s s d r a w b a c k of W h i t t a k e r ' s m e t h o d . T h e r e is, however, a p r o m i s i n g aspect arising f r o m this a p p r o a c h , which is the p o s s i b i l i t y of a p p l y i n g r e c e n t results [32] on a d i a b a t i c j u m p s after separatrix crossing. This has m u c h relevance for the p r o g r a m of f u t u r e studies o u t l i n e d in section 6.

Note added in proof W h i l e this p a p e r was in press, it c a m e to our a t t e n t i o n the w o r k b y Jacques H e n r a r d a n d A n n e L e m a i t r e , " A P e r t u r b a t i v e T r e a t m e n t of the 2 / 1 J o v i a n R e s o n a n c e " . Their analysis i n d i c a t e that W i s d o m ' s effect is s m a l l e r t h a n the eccentricity increase due to t r u n c a t i o n effects. W e observe, however, t h a t the a p p r o a c h w e suggest in section 4.3 differs f r o m t h a t used b y H e n r a r d a n d L e m a i t r e (see their r e m a r k s in the e n d of section 4). W e suspect that our a p p r o a c h m a y e n h a n c e W i s d o m ' s effect.

References [1] a) D. Kirkwood, Meteoric Astronomy (Lippincott, Philadelphia, 1867). b) E.W. Brown, Resonances in the Solar System, Bull. A.M.S. 34 (1928) 265-289. [2] H. Scholl, Resonances in the asteroidal belt, in: Resonances in the Motion of Planets, Satellites and Asteroids, S. Ferraz-Mello and W. Sessin, eds., (Sao Panlo Univ. Press, Brazil, 1985). [3] J. Wisdom, a) The origin of the Kirkwood gaps: a mapping for the asteroidal motion near the 3/1 commensurability, Astr. J. 87 (1982) 577-593. b) Chaotic behavior and the origin of the 3/1 Kirkwood gap, Icarus 56 (1983) 51-74; c) A perturbative treatment of motion near the 3/1 commensurability, Icarus 63 (1985) 272-286. [4] C.D. Murray and K. Fox, Structure of the 3/1 Jovian resonance: a comparison of numerical methods, Icarus 59 (1984) 221-233. [5] H. Poincar6, Les Methodes Nouvelles de la Mecanique Celeste (Gauthier-Villars, Paris, 1897). [6] G.W. Hill, Illustrations of periodic solutions in the problem of three bodies, Astr. J. 22 (1902) 93-97 and 117-121. [7] B.V. Chirikov, A universal instability of many-dimensional oscillator systems, Phys. Rep. 52 (1979) 263-379. [8] B. van tier Pol, On "relaxation-oscillations", Phil. Mag. 2 (1926) 978-992. [9] J. Koiller, J.M. Balthazar and T. Yokoyama, On relaxation chaos: an example from Celestial Mechanics, Proc. First Int. Conf. Physics of Phase-Space, Maryland, 1986, to appear in Springer Lect. Notes in Physics series. [10] a) F. Takens, Implicit differential equations, some open problems, in: Springer Lect. Notes Math. 535, O. Burlet and F. Ronga eds. (Springer, Berlin, 1976). b) P. Cartier, Perturbations singulieres des ~quations differentielles ordinaires et analyse non-standard, Sem. Bourbaki, 1981/1982, n. 580. c) M. Diener, The Canard unchained or how fast/slow dynamical systems bifurcate, Math. Intell. 6:3 (1984) 38-49. d) F. Diener, Sauts des solutions des ~quations E; =f(t, x, 5c), SIAM J. Math. Analysis 17:3 (1986) 533-559. [11] A.T. Sinclair, Periodic solutions close to commensurabilities in the three body problem, Mort. Not. Roy. Astr. Soc. 148 (1970) 325-351.

122

J. Koiller et al./ Relaxation-chaos phenomena in celestial mechanics I

[12] a) V. Arnoid, Mathematical Methods of Classical Mechanics (Springer, Berlin, 1978). b) J.A. Sanders and F. Verhulst, Averaging Methods in Nonlinear Dynamical Systems (Springer, Berlin, 1985). [13] P. Holmes, J. Marsden and J. Scheurle, Exponentially small splitting of separatrices, Center for Pure and Applied Mathematics, Berkeley, PAM-364, 1987. [14] J. Scheurle, A chaotic solution of systems with almost periodic forcing, ZAMP 37:1 (1986). [15] J.A. Sanders, Melnikov's method and averaging, Cel. Mech. 28 (1982) 171-181. [16] J. Koiller, Some simple mechanical systems exhibiting horseshoes and Arnold diffusion, Mat. Aplic. Comp. Brazil 2 (1983) 219-235. [17] H. Aref, Stirring by chaotic advection, J. Fluid Mech. 143 (1984) 1-21. [18] C. Murray, Structure of the 2 : 1 and 3 : 2 Jovian resonances, Icarus 65 (1986) 70-82. [19] R. Leacock and P. O'Connor, Action variable theory and classical frequencies, J. Comp. Phys 62 (1986) 164-179. [20] J. Guckenheimer and P.J. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields (Springer, Berlin, 1983). [21] M. Hirsch, C. Pugh, M. Shub, Invariant Manifolds, Springer Lecture Notes in Math. 583 (Springer, Berlin, 1977). [22] R.C. Robinson, Sustained resonance for a nonlinear system with slowly varying coefficients, SIAM J. Math. Analysis, 14 (1983) 847-860. [23] V.V. Kozlov, Integrability and non-integrability in Hamiltonian mechanics, Russian Math. Surveys 38 (1983) 1-76. [24] P.J. Holmes and J.E. Marsden, Horseshoes in perturbations of Hamiltonian systems with two degrees of freedom, Commun. Math. Phys. 82 (1982) 523-544. [25] P. Veerman and P.J. Holmes, Resonance bands in a two degree of freedom system, Physica 20D (1986) 413-422. [26] P. Byrd and M. Friedman, Handbook of Elliptic Integrals (Springer, Berlin, 1954). [27] a) L. Lerman and Ia. Umanskii, On the existence of separatrix loops in four dimensional systems similar to the integrable Hamiltonian systems, PMM (Applied Math. Mechanics) 47 (1984) 335-340. b) R.C. Robinson, Horseshoes for Hamiltonian systems using the Melnikov integral, Northwestern Univ. Math. Dept. preprint, 1985. [28] Y. Kifer, Ergodic Theory of Random Transformations (Birkhauser, Boston, 1985). [29] R.P. Feynman, Adventures of a Curious Character (Bantam, 1985). [30] S. Kogelman and J.B Keller, Asymptotic theory of nonlinear wave propagation, SIAM J. Applied Math. 24 (1973) 351-352. [31] E.T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies (Princeton Univ. Press, Princeton, 1965). [32] J.L. Tennyson, J.R. Cary and D.F. Esconde, change of the adiabatic invariant due to separatrix crossing, Phys. Rev. Lett. 56 (1986) 2117-2120.