Dynamics in first-order mean motion resonances: analytical study of a simple model with stochastic behaviour
DDynamics in first-order mean motion resonances: analytical study of asimple model with stochastic behaviour
S. Efimov · V. SidorenkoAbstract
We examine a 2DOF Hamiltonian system, which arises in study of first-order mean motionresonance in spatial circular restricted three-body problem “star-planet-asteroid”, and point out somemechanisms of chaos generation. Phase variables of the considered system are subdivided into fast andslow ones: one of the fast variables can be interpreted as resonant angle, while the slow variables areparameters characterizing the shape and orientation of the asteroid’s orbit. Averaging over the fast motionis applied to obtain evolution equations which describe the long-term behavior of the slow variables. Theseequations allowed us to provide a comprehensive classification of the slow variables’ evolution paths. Thebifurcation diagram showing changes in the topological structure of the phase portraits is constructed andbifurcation values of Hamiltonian are calculated. Finally, we study properties of the chaos emerging in thesystem.
Keywords
Hamiltonian system · averaging method · mean motion resonance · chaotic dynamics The model system which will be considered below arises in studies of first-order mean motion resonances(MMR) in restricted three-body problem (R3BP) “star-planet-asteroid”. If asteroid makes p ∈ N revolu-tions around the star in the same amount of time in which the planet makes p + 1 revolutions, there isan exterior resonance of the first-order denoted as p : ( p + 1). The term exterior refers to the fact that inthis case the asteroid’s semi-major axis is larger than semi-major axis of the planet. The interior MMR( p + 1) : p takes place when asteroid makes p + 1 revolutions during the time in which planet makes p . The first-order MMRs are quite common and, therefore, intensively studied by many specialists. Therelated bibliography is given in Gallardo (2018) and Nesvorny et al. (2002). In particular, much effort hasbeen spent to reveal why 2 : 1 resonance with Jupiter corresponds to one of the largest gaps in the mainasteroid belt (so-called Hecuba gaps), whereas 3 : 2 MMR resonance is populated by numerous objectsof Hilda group, and it is also very likely, that Thule group of objects in 4 : 3 MMR is rather large (Broz S. EfimovMoscow Institute of Physics and Technology9 Institutskiy per., Dolgoprudny, Moscow Region, 141701, Russian FederationE-mail: efi[email protected]. SidorenkoKeldysh Institute of Applied MathematicsRussian Academy of Sciences,Miusskaya Sq., 4, 125047 Moscow, Russian FederationandMoscow Institute of Physics and Technology9 Institutskiy per., Dolgoprudny, Moscow Region, 141701, Russian Federation a r X i v : . [ phy s i c s . s p ace - ph ] D ec and Vokrouhlicky 2008; Henrard 1996; Lemaitre and Henrard 1990). The discoveries of trans-Neptunianobjects made it urgent to study the exterior resonances with Neptune: twotino (MMR 1 : 2) and plutino (MMR 2 : 3) form big subpopulations in the Kuiper belt (Li et al. 2014a,b; Nesvorny and Roig 2000,2001).It is possible to construct a model of dynamics in first-order MMR, taking into account only the leadingterms in the Fourier series expansion of disturbing function (Sessin and Ferraz-Mello 1984; Wisdom 1986;Gerasimov and Mushailov 1990). However, studies of planar R3BP (Beaug´e 1994; Jancart et al. 2002)revealed, that some important characteristics of first-order MMRs are reproduced only when the second-order Fourier terms are accounted for. In this paper we concentrate our attention on that part of a phasespace, where eccentricities and inclinations are in relation e (cid:28) i (cid:28) e > i = 0). We intend to demonstrate, that innon-planar case second-order terms are no less important, as they make a model essentially stochastic.In contrast, the first-order models are proven to be integrable (Sessin and Ferraz-Mello 1984), and thuscannot reproduce chaotic dynamics found in multiple numerical studies of first-order resonances (Wisdomand Sussman 1988; Giffen 1973; Winter and Murray 1997a,b; Wisdom 1987),There are different mechanisms for generating chaos in the dynamics of celestial bodies (Holmes 1990;Lissauer 1999; Morbidelli 2002). Presence of MMR may lead to the so-called adiabatic chaos (Wisdom1985), which is caused, roughly speaking, by small quasi-random jumps between regular phase trajectoriesin certain parts of the phase space, where adiabatic approximation is violated. Applying systematicallyWisdom’s ideas to study of MMRs (Sidorenko 2006; Sidorenko et al. 2014; Sidorenko 2018), we foundthat adiabatic chaos often coexists with the quasi-probabilistic transitions between specific phase regions.Both phenomena occur in that part of the phase space, where the “pendulum” or first-order SecondFundamental Model for Resonance (Henrard and Lamaitre 1983) approximations fail, because the firstharmonic in the disturbing function Fourier series is not dominant. The goal of this paper is to carry outa comprehensive analysis of the introduced second-order model and investigate described mechanisms ofchaotization, which, in our opinion, have not received proper attention in the past.The paper is organized as follows. In Section 2 the model Hamiltonian system, which has a structureof slow-fast system, is introduced. In Section 3 the fast subsystem is studied. Equations of motion for slowsubsystem are constructed in Section 4 and their solutions are analyzed in Section 5. Section 6 is devotedto different chaotic effects present in the discussed model. In Section 7 numerical evidence for existence ofdescribed phenomena is shown. The results are summarized in the last section. In Appendix A we revealhow the proposed model was derived. Details of the averaging procedure are elucidated in Appendix B. We are dealing with 2DOF Hamiltonian systems with specific symplectic structure: dϕdτ = ∂Ξ∂Φ , dΦdτ = − ∂Ξ∂ϕ ,dxdτ = ε ∂Ξ∂y , dydτ = − ε ∂Ξ∂x . (1)The Hamiltonian Ξ in (1) is expressed by Ξ ( x, y, ϕ, Φ ) = Φ W ( x, y, ϕ ) , W ( x, y, ϕ ) = x cos ϕ + y sin ϕ + cos 2 ϕ. (2)Appendix A describes in detail how the system (1)-(2) arises in studies of first-order MMR in three-dimensional R3BP. Here we only note that ε ∼ µ / , x ∝ e cos ω, y ∝ − e sin ω, where µ (cid:28) e and ω denotethe eccentricity and the argument of pericenter of asteroid’s osculating orbit respectively. Further ε istreated as a small parameter of the problem. Since in general variables ϕ, Φ, x, y vary with different rates ( dϕ/dτ, dΦ/dτ ∼
1, while dx/dτ, dy/dτ ∼ ε (cid:28) fast subsystem (described bythe first line of equations) and slow subsystem (the second line).In limiting case ε = 0 equations of fast subsystem coincide with equations of motion for particle withunit mass in a field with potential W ( x, y, ϕ ), where x and y are treated as parameters. Let ϕ ( τ, x, y, ξ ) , Φ ( τ, x, y, ξ ) (3)be a solution of fast subsystem with fixed values of x , y , and value ξ of Hamiltonian Ξ , which is the firstintegral of system (1). In general the resonant angle ϕ in (3) can librate between two constant values orchange monotonously through the whole interval [0 , π ), i.e. circulate . In either case ϕ ( τ + T, x, y, ξ ) = ϕ ( τ, x, y, ξ ) mod (2 π ) , where T ( x, y, ξ ) is a period of the solution (3).Because fast variables vary much faster then the slow ones, the right-hand sides in differential equationsof slow subsystem in (1) can be replaced by their average values along the solution (3). This yields the evolution equations , which describe secular variations of x and y in closed form: dxdτ = ε (cid:28) ∂Ξ∂y (cid:29) , dydτ = − ε (cid:28) ∂Ξ∂x (cid:29) . (4)Here (cid:104) Λ (cid:105) = 1 T ( x, y, ξ ) T ( x,y,ξ ) (cid:90) Λ ( x, y, ϕ ( τ, x, y, ξ )) dτ . (5)The solution (3) has an action integral: J ( x, y, ξ ) = 12 π T ( x,y,ξ ) (cid:90) Φ ( τ, x, y, ξ ) dτ . (6)For ε (cid:54) = 0 function J ( x, y, ξ ) becomes an adiabatic invariant of slow-fast system (1). For a fixed ξ trajec-tories of averaged equations (4) on a phase plane ( x, y ) go along the lines with constant values of J . Thisallows to classify evolution equations (4) as adiabatic approximation (Neishtadt 1987a; Wisdom 1985).In the next Sections we go through all the steps in construction of evolution equations via describedapproach and analyze in detail the behaviour of slow variables on different levels of Hamiltonian Ξ = ξ . ε = 0) x, y ) based on the number of librating solutionsBecause the variables x and y change very slowly (1), they can be treated as constant parameters, whenconsidering the motion of the fast subsystem. Then the potential W is just a “two-harmonic” function of ϕ defined on circle S , and the motion in such potential can be described in terms of elliptic functions.There are different types of motion depending on the Hamiltonian level ξ at which it occurs (Fig. 1).For us it is important, that for some values of x and y two different librating solutions can exist on thesame ξ level (Fig. 1c). This situation can take place when W ( ϕ ) has four extrema on S .A necessary condition for extremum ∂W / ∂ϕ = 0 after the replacement λ = tan ( ϕ/
2) yields yλ + 2( x + 4) λ + 2( x − λ − y = 0 . (7)Let A denote a region on the plane ( x, y ), in which W ( ϕ ) has four extrema. The equation (7) has four realroots inside this region and only two outside. Thus on the border of the region A the number of unique Fig. 1
Levels ξ of Hamiltonian Ξ corresponding to different types of fast subsystem’s motion: a. circulation, b.libration, c. two coexisting librating solutions, d. the motion is impossible Fig. 2
Extremal surface of potential W ( x, y, ϕ ) (left), and astroid bounding the region A , where W has fourextrema as a function of ϕ on S (right). A , . . . , A – astroid’s cusps real roots is 3 (with the exception of finite number of points in which there is only one unique real root)and the discriminant of (7) is equal to zero. Therefore the equation for the border of the region A is x + 3 x y − x + 3 x y + 336 x y + 768 x + y − y + 768 y − . (8)By collecting the parts of this equation into perfect cube (8) is transformed to canonical algebraic equationof astroid (Fig. 2): (cid:16) x + y − (cid:17) + 27 · x y = 0 . (9)Which can be further reduced to x / + y / = 4 / . (10) It is convenient to use this astroid for the reference on the phase plane ( x, y ).3.2 Critical curve partitioning the plane ( x, y ) into regions with different types of fast subsystem’ motionLet us introduce the notations W min ( x, y ) and W max ( x, y ) for global minimum and maximum of function W on S for given values of slow variables. If ( x, y ) ∈ A , then W has the second pair of minimum andmaximum, which we shall denote W ∗ min ( x, y ) and W ∗ max ( x, y ) respectively. Using these auxiliary functions,we can partition the ( x, y ) plane for a given ξ into different regions based on the type of fast subsystem’smotion: Q = { ( x, y ) | ξ < W min } ,Q = (cid:8) ( x, y ) / ∈ A (cid:12)(cid:12) ξ ∈ (cid:0) W min , W max (cid:1) (cid:9) (cid:91)(cid:8) ( x, y ) ∈ A (cid:12)(cid:12) ξ ∈ (cid:0) W min , W max (cid:1) \ (cid:0) W ∗ min , W ∗ max (cid:1) (cid:9) ,Q = (cid:8) ( x, y ) ∈ A (cid:12)(cid:12) ξ ∈ (cid:0) W ∗ min , W ∗ max (cid:1) (cid:9) ,Q = { ( x, y ) | ξ > W max } . The region Q ( ξ ) will be called a forbidden region , because inside of it Ξ < ξ for any values of fast variablesand fast subsystem has no solutions (Fig. 1d). Region Q ( ξ ) is the region with the single librating solution(Fig. 1b), Q ( ξ ) is the region with two librating solutions at given level ξ (Fig. 1c), and finally the region Q ( ξ ) is where fast subsystem’s solution circulates (Fig. 1a). Illustrations for regions Q ( ξ ) , . . . , Q ( ξ ) willfollow. Fig. 3
Tangency of Hamiltonian level ξ and different extrema of W ( ϕ ), which occurs on the borders of regions Q ( ξ ) , . . . , Q ( ξ ) Before that let us consider a border Γ ( ξ ) between these regions. In every point of the border value of W is equal to ξ in one of its critical points (cf. Figures 1 and 3), which is why we shall call Γ ( ξ ) a criticalcurve . After replacement λ = tan ( ϕ/
2) the equation W ( ϕ ) = ξ transforms into algebraic equation:(1 − ξ − x ) λ + 2 yλ − ξ + 3) λ + 2 yλ + ( x + 1 − ξ ) = 0 . (11) As Figure 3 demonstrates, the point ( x, y ) lies on the critical curve when equation (11) have at least onemultiple real root, which is equivalent to discriminant of (11) being equal to zero: D ( x, y, ξ ) = 64 ξ − ξ − x + ξ x − ξx − x y − x ++ 16 ξ x − ξ x − ξx − x y + 2 ξ x y + 78 x y − x −− y + ξ y + 18 ξy − y − ξ y − ξ y + 144 ξy − y + 64 = 0 . (12)Thus in the regions Q ( ξ ), Q ( ξ ), Q ( ξ ) the discriminant D ( x, y, ξ ) >
0, while D ( x, y, ξ ) < Q ( ξ ),and D ( x, y, ξ ) = 0 on the critical curve Γ ( ξ ). Figure 4 depicts Γ ( ξ ) on the plane of slow variables fordifferent values of ξ . At | ξ | < − < ξ < −
1, the curve intersects itself on axis x , with the x coordinates of self-intersection points being defined by equation x + 8( ξ + 1) = 0 . If 1 < ξ <
3, points of self-intersection lie on y axis, and their y coordinates are defined by equation y − ξ −
1) = 0 . Fig. 4
Critical curve Γ ( ξ ): shape of the curve for different ξ values (left), and parametrization of the curve by (cid:95) ϕ with arrows showing the direction in which the parameter increases (right) The critical curve allows a parametrization Γ ( ξ ) = (cid:110) x = cos (cid:95) ϕ ( ξ + cos 2 (cid:95) ϕ − , y = sin (cid:95) ϕ ( ξ + cos 2 (cid:95) ϕ + 2) (cid:12)(cid:12) (cid:12) (cid:95) ϕ ∈ S (cid:111) , (13)which is illustrated by Figure 4. The parameter (cid:95) ϕ in (13) coincide with the critical points of potential W ( ϕ ) at given level ξ , as depicted in Figure 3, in respective points of ( x, y ) plane: W ( (cid:95) ϕ, x ( (cid:95) ϕ, ξ ) , y ( (cid:95) ϕ, ξ )) = ξ, ∂∂ϕ W ( ϕ, x ( (cid:95) ϕ, ξ ) , y ( (cid:95) ϕ, ξ )) (cid:12)(cid:12)(cid:12) ϕ = (cid:95) ϕ = 0 . The introduced parametric representation is convenient, in particular, for defining the location ofcusps and self-intersection points. For cusps (cid:95) ϕ = ϕ ∗ , where ϕ ∗ is obtained from the equation tan ϕ ∗ =(3 + ξ )/(3 − ξ ), which have four roots on S . We shall denote these cusps as Y ,..., Y with the lowerindex being the number of a quadrant, in which the respective value ϕ ∗ lies. Self-intersection points of Γ ( ξ ) on x axis ( − < ξ < −
1) we shall denote as B and B for right and left half-planes respectively.Self-intersection points on y axis (1 < ξ <
3) we shall denote as S and S for upper and lower half-planes. Note:
It can be demonstrated, that curves Γ ( ξ ) are the involutes of astroid (10) constructed withtethers of length 3 ± ξ extended from astroid’s cusps. This makes Γ ( ξ ) also a family of equidistant curveswith the distance | ξ a − ξ b | between any two curves Γ ( ξ a ) and Γ ( ξ b ).3.3 Transformation of regions Q ( ξ ) , ..., Q ( ξ ) with change of ξ It should be noted first, that region with a single librating solution Q ( ξ ) is present on plane ( x, y ) for allvalues of ξ . Other regions appear and disappear, as ξ crosses several bifurcation values ξ i : ξ = − , ξ = − , ξ = 1 , ξ = 3 . We shall describe, how the regions are transformed, as ξ increases.If ξ < ξ , there exists a forbidden region Q ( ξ ) around the point (0 , x, y ) planebeing the Q ( ξ ) region.At ξ = ξ on the right and on the left from region Q ( ξ ) two parts of region Q ( ξ ) (region with twolibrating solutions) appear (Fig. 5).At ξ = ξ the region Q ( ξ ) disappears and Q ( ξ ) becomes connected (Fig. 6).At ξ = ξ region Q ( ξ ) is separated into two parts again by appearing region Q ( ξ ) (region withcirculating resonant angle) around the point (0 ,
0) as seen in Figure 7.At ξ = ξ region Q ( ξ ) disappears (Fig. 8). For ξ > ξ there exists only region Q ( ξ ) surrounded by Q ( ξ ). Fig. 5
Bifurcation at ξ = ξ : appearance of region Q . Here and further region Q is colored dark gray, region Q – orange. The rest blank space corresponds to region Q Fig. 6
Bifurcation at ξ = ξ : vanishing of forbidden region Q Fig. 7
Bifurcation at ξ = ξ : appearance of region Q (here and further colored green) Fig. 8
Bifurcation at ξ = ξ : vanishing of region Q Let us now describe how borders of regions Q ( ξ ) , ..., Q ( ξ ) transforms with increase of ξ . The borderbetween Q ( ξ ) and Q ( ξ ) we shall call the existence curve and denote it as Γ ( ξ ). It corresponds to thepart of the critical curve Γ ( ε ) in which W min ( x, y ) = ξ . For ξ < ξ curve Γ ( ε ) ≡ Γ ( ε ). For ξ < ξ < ξ curve Γ ( ε ) consists of two intervals of Γ ( ε ), lying between points of self-intersection B and B .We shall adopt the traditional terminology common in studies of slow-fast systems (Wisdom 1985;Neishtadt 1987a) with modifications made to better represent the specifics of the discussed problem.The border between regions Q ( ξ ) and Q ( ξ ) we shall call an uncertainty curve of the first kind anduse a notation Γ ( ξ ) for it. Points of the uncertainty curve of the first kind are defined by condition W max ( x, y ) = ξ . If ξ > ξ , then Γ ( ξ ) ≡ Γ ( ξ ). For ξ < ξ < ξ , the curve Γ ( ξ ) consists of Γ ( ξ ) parts,which are contained between points of self-intersection S and S .The part of the border between Q ( ξ ) and Q ( ξ ), along which holds the equality W ∗ max ( x, y ) = ξ ,we shall call an uncertainty curve of the second kind and denote it Γ ( ξ ). For ξ < ξ < ξ the Γ ( ξ ) = Y Y ∪ Y Y , where Y Y and Y Y are segments of Γ ( ξ ), which lie between corresponding cusps. If ξ < ξ < ξ then curve Γ ( ξ ) = S Y ∪ S Y ∪ S Y ∪ S Y . For the rest part of the border between Q ( ξ )and Q ( ξ ) holds W ∗ min ( x, y ) = ξ . As no dynamical effects of interest are happening on this segment, weshall not refer to it further.Figure 9 depicts a diagram, that shows values of (cid:95) ϕ defining positions of cusps and self-intersectionpoints on Γ ( ξ ), as well as the segments which correspond to existence curve and two uncertainty curves.Due to the symmetry, it is sufficient to consider only (cid:95) ϕ ∈ [0 , π/
2] (Fig. 4).3.4 Three-dimensional representation of the set of curves Γ ( ξ )Curves Γ ( ξ ) can be interpreted as cross sections of some surface F in the space xyξ by equi-Hamiltonianplanes ξ = const (Fig. 10a,b). In this space for fixed value of (cid:95) ϕ the equations (13) define a straight line,which means that the surface F is ruled (Fig. 10c). Fig. 9
Diagram showing the partition of critical curve into existence curve Γ ( ξ ) and critical curves Γ , ( ξ ) Fig. 10
Surface F composed of curves Γ ( ξ ) in xyξ space: a. general representation of the surface, b. surface’s crosssections by equi-Hamiltonian planes, c. rulings of surface F The same surface F defined by the equation analogous to (12) also appears in a completely differentproblem studied by Batkhin (2012), where it partitions the parametric space of some mechanical systeminto regions with different stability properties. ∂Ξ∂x = cos ϕ, ∂Ξ∂y = sin ϕ, construction of the evolution equations (4) require calculating two averaged properties: (cid:104) sin ϕ (cid:105) = T ( x,y,ξ ) T ( x,y,ξ ) (cid:82) sin ϕ ( τ, x, y, ξ ) dτ, (cid:104) cos ϕ (cid:105) = T ( x,y,ξ ) T ( x,y,ξ ) (cid:82) cos ϕ ( τ, x, y, ξ ) dτ. (14)The equality dϕdτ = Φ = ± (cid:112) ξ − W ( x, y, ϕ )] allows finding a period of fast subsystem’s solution at Hamiltonian level Ξ = ξ and calculating (afterproper change of variables) values of integrals on the right-hand side of (14), e.g. for librating solutions: T ( x, y, ξ ) = 2 ϕ ∗ (cid:82) ϕ ∗ dϕ √ ξ − W ( x,y,ϕ )] , T ( x,y,ξ ) (cid:82) f ( ϕ ( t, x, y, ξ )) dτ = 2 ϕ ∗ (cid:82) ϕ ∗ f ( ϕ ) dϕ √ ξ − W ( x,y,ϕ )] . (15)Here ϕ ∗ and ϕ ∗ denote minimum and maximum values of angle ϕ in librating solution, along whichthe averaging is being performed.For the system (2) period T ( x, y, ξ ) and integrals (14) can be expressed in terms of complete ellipticalintegrals of the first and the third kind. Further the concise description of this transformation is given,using the case − π ≤ ϕ ∗ < ϕ ∗ ≤ π (16)as an example. After standard substitution λ = tan ( ϕ/
2) we obtain T ( x, y, ξ ) = 4 λ ∗ (cid:82) λ ∗ dλ √ R ( λ ) , (cid:104) sin ϕ (cid:105) = ϕ ∗ (cid:82) ϕ ∗ sin ϕdϕ √ ξ − W ( x,y,ϕ )] =4 λ ∗ (cid:82) λ ∗ λdλ (1+ λ ) √ R ( λ ) , (cid:104) cos ϕ (cid:105) = ϕ ∗ (cid:82) ϕ ∗ cos ϕdϕ √ ξ − W ( x,y,ϕ )] =2 λ ∗ (cid:82) λ ∗ ( − λ ) dλ (1+ λ ) √ R ( λ ) . (17)Here λ ∗ = tan ( ϕ ∗ / λ ∗ = tan ( ϕ ∗ / R ( λ ) in (17) is a quartic polynomial R ( λ ) = d λ + d λ + d λ + d λ + d with coefficients d = ξ − x, d = d = − y, d = 2 ξ + 6 , d = ξ − − x. Integrals on the right-hand side in (17) can be rewritten as follows: T ( x, y, ξ ) = √ | d | I , , ϕ ∗ (cid:82) ϕ ∗ sin ϕdϕ √ ξ − W ( x,y,ϕ )] = √ | d | I , , ϕ ∗ (cid:82) ϕ ∗ cos ϕdϕ √ ξ − W ( x,y,ϕ )] = √ | d | (cid:0) I , − I , (cid:1) , (18)where notation I k,r is used for integrals: I k,r = λ ∗ (cid:90) λ ∗ λ r dλ ( λ + 1) k (cid:112) ± ( λ − a )( λ − a )( λ − a )( λ − a ) . (19) a k – roots of polinomial R ( λ ), and the sign in the square root is determined by sign of coefficient d .Integrals (19) can be presented as linear combinations of elliptic integrals of the first and the third kind: I , = c , K ( k ) ,I , = c , K ( k ) + c , Π ( h, k ) + ¯ c , Π (cid:0) ¯ h, k (cid:1) ,I , = g , K ( k ) + g , Π ( h, k ) + ¯ g , Π (cid:0) ¯ h, k (cid:1) . (20)Analytical expressions for coefficients c m,l , g m,l , moduli k , and parameters h in (20) depend on integrationinterval and properties of R ( λ ) roots. These expressions are gathered in the Appendix B. Expressions of the kind cΠ ( h, k ) + ¯ cΠ (¯ h, k ) ( c, h ∈ C , k ∈ R , < k < λ = tan ( ϕ/
2) the integration in (17) is carried over two semi-infinite intervals. E.g., for ϕ ∗ < π, ϕ ∗ > π, ϕ ∗ − ϕ ∗ < π the expression for the period T is T ( x, y, ξ ) = 4 λ ∗ (cid:90) −∞ dλ (cid:112) R ( λ ) + + ∞ (cid:90) λ ∗ dλ (cid:112) R ( λ ) . In these cases it is implied that all I k,r in (18) are the sum of two integrals as well.After all described transformations evolution equations (4) take a simple form dxdτ = ε I , I , , dydτ = ε (cid:20) − I , I , (cid:21) . (21)It should be noted, that there is an ambiguity in calculation of the right-hand side parts of evolutionequations in the region Q ( ξ ): the result depends on the choice of the fast subsystem’s solution, and thereare two different librating solutions in Q ( ξ ). Consequently, phase portraits of (4) shall contain two setsof trajectories in Q ( ξ ), which correspond to two possible variants of slow variables’ evolution.When describing the crossing of uncertainty curves Γ i ( ξ ) by the projection ζ ( τ ) = ( x ( τ ) , y ( τ )) T ofthe phase point z ( τ ) = ( ϕ ( τ ) , Φ ( τ ) , x ( τ ) , y ( τ )) T on the plane ( x, y ), we shall confine ourselves to formalcontinuation of averaged system’s trajectories, lying on opposite sides of Γ i . The detailed analysis of theseevents is given in Neishtadt (1987a), Neishtadt and Sidorenko (2004), and Sidorenko et al. (2014).4.2 Fast subsystem’s action variable – integral of evolution equationsAs was mentioned in the Section 2, adiabatic invariant J ( x, y, ξ ) of slow-fast system is a first integral ofevolution equations (4). Formula (6) for J ( x, y, ξ ) can be expressed as a linear combination of integrals I k,r defined by (19): J ( x, y, ξ ) = 1 π ϕ ∗ (cid:90) ϕ ∗ Φdϕ = 1 π ϕ ∗ (cid:90) ϕ ∗ (cid:112) ξ − W ( x, y, ϕ )] dϕ = 2 π λ ∗ (cid:90) λ ∗ (cid:112) R ( λ ) dλ (1 + λ ) == 2 √ π (cid:112) | d | (cid:2) d I , + d I , + ( d − d ) I , + ( d + d − d ) I , (cid:3) . (22)Analytical expression for I , is presented in the Appendix B alongside the rest of the integrals previouslyshown in (20). Fig. 11
Phase portraits of evolution equations:a. ξ = −
4, b. ξ = − .
5, c. ξ = − .
1, d. ξ = 1 .
45, e. ξ = 3 .
4, f. ξ = 8 . ξ < ξ the structureof phase portrait is simple – all phase trajectories are represented by closed loops encircling the forbiddenregion Q ( ξ ) (Fig. 11a). Figure 11b depicts a typical phase portrait for ξ ∈ ( ξ , ξ ) – two symmetrical partsof region Q ( ξ ) adjoining the central region Q ( ξ ) contain two layers of phase trajectories. For ξ ∈ ( ξ , ξ )there are only two regions on the phase plane, i.e. Q ( ξ ) and Q ( ξ ) (Fig. 11c). In addition to presentedin Figure 11c,d change of phase portrait’s global structure, the behaviour of phase trajectories near theuncertainty curves Γ , ( ξ ) at ξ ∈ ( ξ , ξ ) have some specific qualitative differences at different ξ values.The detailed description of that is given in the end of this section.Phase portraits at ξ > ξ have five stationary points: the origin of xy -plane is the stable point ofthe center type, two more center points are symmetrically located on the y -axis, and two unstable saddlepoints are symmetrically located on the x -axis (Fig. 11d–f). Phase portraits depicted in Figures 11e and11f are differ in relative positions of heteroclinic trajectories and uncertainty curve Γ .Ordinates y of the center points above and below plane’s origin are defined by equation K ( m ) = 2 Π ( n | m ) , (23)where m = U + U − , n = U + ξ + 1 − | y | , U ± = ξ − ± (cid:112) y + 8(1 − ξ ) . Abscissae x of the saddle points are defined by the same relation (23) with different parameters: m = Q + − Q −− Q ++ Q − + , n = Q + − Q −− , Q ±± = (cid:112) x + 8(1 + ξ ) ± ± | x | . Solutions of (23) are plotted in Figure 12 as functions of ξ for both types of stationary points. Note,that top and bottom center points are located in Q at ξ < ξ = 2 and in Q at ξ > ξ (Fig. 13). Fig. 12
Coordinates of evolution equations’ stationary points, which lie on axis x (left), and on axis y (right).Values of coordinates as functions of ξ are represented by violet lines. The rest of the coloring is consistent withFigure 11 in denoting the regions Q i and points from different parts of the critical curve4 Fig. 13
Transition of the top center point from Q ( ξ ) to Q ( ξ ) at ξ = ξ = 2 ξ value, points of uncertaintycurves’ intersections with limiting trajectories, that are limiting cases for different families of trajectories,should be considered. We shall refer to them as limiting points . Fig. 14
Limiting points R i and I i at ξ ≈ Several types of limiting points are depicted in Figure 14. Trajectories, which go in Q between twopoints on uncertainty curve Γ , can be divided into two families: the trajectories that intersect x axis, andthose that do not. Thus there is limiting trajectory that separates these two families (in Figure 14 it iscolored red). We shall denote the limiting points corresponding to this trajectory as R , ..., R (Fig. 14). The subset of trajectories, which do not cross x -axis in Q as another limiting case contain trajectories,which come to uncertainty curve from the side of Q and then reflect back without exiting to Q . InFigure 14 these trajectories are colored purple, and the corresponding limiting points are denoted I , ..., I . Fig. 15
Limiting points K i , M i and V i ( ξ ∼ . More limiting points are depicted in figure 15. Points K , ..., K ∈ Γ are connected to cusps Y , ..., Y by limiting phase trajectories, which separate the family of trajectories lying to the one side of uncertaintycurve Γ from those which in Q connect symmetric points on left and right sides of uncertainty curves Γ , . Limiting points M , ..., M divides each of four Γ segments (located in each of four quadrants of xy plane) into two sections: one section has both trajectories, which adjoin from Q side, going in thesame direction, while trajectories adjoining the other section go in opposite directions (also see Figure 19).Points V , ..., V of Γ curve’s intersections with mentioned earlier separatrices constitute the last type oflimiting points.All differences between the phase portraits (Fig. 11) emerge from changes in position of limiting pointson the critical curve. By adding these points to diagram in Figure 9 we obtain a bifurcation diagram(Fig. 16), from which all changes in topological structure of phase portraits can be understood. Let usdescribe, what is happening to different limiting points by going successively from low to high values of ξ . The first limiting points – R i and I i – appear at ξ = ξ ≈ − . ξ = ξ ≈ . R i merge with cusps Y i , and at the same time points K i appear. Points M i and V i emerge simultaneouslywith region Q and self-intersections points S i of the critical curve at ξ = ξ . All points on the Γ part of the critical curve – I i , K i , and M i , as well as the ends Y i , S i of Γ itself – disappear at ξ = ξ mergingwith the astroid’s cusps A and A . Finally, at ξ = ξ ≈ . V i disappear by merging pairwise. Fig. 16
Bifurcation diagram showing the dependance of limiting points position on the critical curve Γ on ξ . Bifurcation Hamiltonian values ξ , . . . , ξ partition the whole range of ξ into nine intervals, meaningthat there is a total number of nine different types of phase portraits for slow subsystem. The bifurcationat ξ (Fig. 13) is omitted from Figure 16 because it does not bear any significance for the dynamical effectsconsidered further. Γ ( ξ ) three guiding trajectories , i.e. trajectories of averaged system (4), meet: two onthe side of Q ( ξ ) region and one on the Q ( ξ ) side. In the case, when two out of these three trajectoriesare outgoing, the transition of the phase point to either one of them can be considered as a proba-bilistic event. In the original system (1) initial values of fast variables corresponding to two differentoutcomes are strongly mixed in the phase space. Therefore even small variation of initial conditions z (0) = ( ϕ (0) , Φ (0) , x (0) , y (0)) T can lead to qualitative change in system’s evolution. As an example, Fig-ure 17 depicts projections on the plane ( x, y ) of two trajectories γ , , obtained as solutions of the system(1) with close initial conditions. Both trajectories approach uncertainty curve along the same guidingtrajectory, but diverge after that – γ exits Q ( ξ ) and follows the outgoing guiding trajectory in Q ( ξ ),while γ turns and goes back along the other guiding trajectory in Q ( ξ ). Fig. 17
Two phase trajectories starting from very close initial conditions may diverge at uncertainty curve Γ , as γ and γ do. Green and red lines show two guiding trajectories departing from the same point of Γ , to which theblue one arrives. Trajectory γ of non-averaged system may go along either one of them after reaching the borderbetween Q and Q In deterministic systems with strongly entangled trajectories the probability of a specific outcome isdetermined by the fraction of phase volume occupied by corresponding initial conditions (formal definitioncan be found in Arnold (1963) and Neishtadt (1987b)). In order to find the probabilities of transitions todifferent outgoing trajectories in some point ( x ∗ , y ∗ ) on Γ , two auxiliary parameters must be calculatedfirst (Neishtadt 1987b; Artemyev et al. 2013): Θ , = + ∞ (cid:90) −∞ (cid:18) ∂W ∗ max ∂x ∂W∂y − ∂W ∗ max ∂y ∂W∂x (cid:19) ϕ s , ( t,x ∗ ,y ∗ ,ξ ) dτ . (24)These parameters have a meaning of rates with which areas bounded by separatrices ( ϕ s ( τ, x ∗ , y ∗ , ξ ) , Φ s ( τ, x ∗ , y ∗ , ξ ))and ( ϕ s ( τ, x ∗ , y ∗ , ξ ) , Φ s ( τ, x ∗ , y ∗ , ξ )) in fast variables’ phase space change. After substitution of specificpotential function (2) into (24) and change of integration variable to ϕ we obtain: Θ = 2 (cid:95) ϕ (cid:90) ϕ s min sin( ϕ − (cid:95) ϕ ) (cid:112) ξ − W ( x ∗ , y ∗ , ϕ )) dϕ, Θ = 2 ϕ s max (cid:90) (cid:95) ϕ sin( ϕ − (cid:95) ϕ ) (cid:112) ξ − W ( x ∗ , y ∗ , ϕ )) dϕ. (25)Here (cid:95) ϕ is the coordinate ϕ of the saddle point in the fast subsystem’s phase portrait (it coincides withthe value (cid:95) ϕ corresponding to point ( x ∗ , y ∗ ) in (13)); ϕ s min and ϕ s max are the minimal and the maximalvalues of ϕ in homoclinic trajectories to the left and to the right of the saddle point respectively. Applyingthe substitution λ = cot[( ϕ − (cid:95) ϕ ) /
2] to (25) we obtain: Θ , = 4 √ A λ max (cid:90) λ min λ dλ (1 + λ ) (cid:112) ( λ − a )( λ − b ) , (26)where a = B − √ B − ACA , b = B + √ B − ACA ,A = ξ + 3 cos(2 (cid:95) ϕ ) , B = 2 sin(2 (cid:95) ϕ ) , C = ξ − cos(2 (cid:95) ϕ ) . It should be noted, that inequalities
A > B > AC hold for all points of Γ . Integration in(26) is carried out over the intervals ( −∞ , a ) and ( b, + ∞ ) for Θ and Θ respectively. The result can beobtained by using Cauchy’s residue theorem: Θ = 4 √ A Re (cid:18) L − iπαβ (cid:19) , Θ = 4 √ A Re (cid:18) Lαβ (cid:19) . Here L = ln (cid:18) b − a ( α − β ) (cid:19) , α = √ a + i, β = √ b + i, and the branches of multifunctions are selected in such way, that Im( L ), arg( α ), and arg( β ) ∈ (0 , π ).Now we can write down the expressions for probabilities of different evolution scenarios for a phasepoint on Γ ( ξ ). Let us denote the probability of point going to region Q ( ξ ) as P , and probabilitiescorresponding to two trajectories going inside Q ( ξ ) as P and P . The resulting formulae (Artemyevet al. 2013) take a form: P = 1 − P − P , P = max( Θ , /Θ Σ , P = max( Θ , /Θ Σ , (27)where Θ Σ = max( Θ ,
0) + max( Θ ,
0) + max( − Θ − Θ , . In Figure 18 and Figure 19 change of probabilities (27) along Γ is shown alongside with phaseportraits at corresponding values of Hamiltonian and limiting points, which mark the change of sign in Θ , or Θ + Θ . Fig. 18
Change of P i along Y Y segment of Γ at ξ = 0. The graph for Y Y can be reconstructed by thesymmetry, changing the y sign and swapping red and blue plots. In I , . . . , I the sum Θ + Θ = 0, which resultsin P plot sticking to 0 on one side from these points Fig. 19
Change of P i along segments Y S and S Y of Γ at ξ = 2 .
4. In M , parameter Θ = 0, and in M , parameter Θ = 0, which results in singularities of probabilities plots in these limiting points9 Adiabatic chaos emerges due to non-applicability of adiabatic approximation near the uncertainty curve.As a result the projection of phase point ζ ( τ ) = ( x ( τ ) , y ( τ )) T leaves the vicinity of uncertainty curve alongthe guiding trajectory, which slightly differs from the direct continuation of approach trajectory (Fig. 20).The resulting offset between incoming and outgoing guiding trajectories can be treated as a quasi-randomjump with order of magnitude ε | ln ε | (Tennyson et al. 1986; Neishtadt 1987b,a).As a result of persistent jumps any trajectory that crosses the uncertainty curve after a long timewill fill a whole region of phase plane, which we will call an adiabatic chaos region (Fig. 21). This regionconsists of points belonging to trajectories, which cross the uncertainty curve. Fig. 20
The jump of trajectory γ from guiding trajectory γ to γ upon crossing the uncertainty curve Fig. 21
Adiabatic chaos region at ξ = 1 . At ξ ∈ ( ξ , ξ ) the uncertainty curve is crossed by the separatrices, which connect two saddle points. Aphase point projection moving along a guiding trajectory close to separatrix, when crossing the uncertainty curve, may jump over the separatrix and begin to move along the other guiding trajectory belongingto completely different family. Thus the properties of long-term evolution are suddenly changed on aqualitative level. E.g., in Figure 21 the motion of a phase point projection ζ ( τ ) circling around thecoordinate origin in the central part of adiabatic chaos region by crossing the separatrix may transforminto circulation in opposite direction around one of two center points (23), which lie on y -axis in upperor lower half-plane. This event can be interpreted as a capture into Kozai-Lidov resonance and it isaccompanied by decrease of average inclination value, about which the long-term oscillations occur. Construction of the discussed analytical model involved several assumptions that may seem loose. It isthus required to test whether the model can be applied to orbital dynamics of real life objects and theseassumptions were not overly restrictive.For this purpose we used Mercury integrator (Chambers 1999) and carried out several numericalsimulations of the Solar system composed of the Sun, the four giant planets, and about 700 known Kuiperbelt objects (KBO) near 1 : 2, 2 : 3, and 3 : 4 resonances with Neptune represented by test particles(masses of four inner planets were added to the Sun in order to facilitate the integration). Total time ofintegration 15
Myr is one order of magnitude larger then the characteristic time T N /ε ≈ . Myr ofslow variables evolution (given the orbital period of Neptune T N ≈ y and Neptune/Sun mass ratiodefining the small parameter ε ∼ − ). Fig. 22
Comparison of the Solar system’s numerical integration with analytical model. Several Kuiper belt objectswith Hamiltonian values close to ξ ≈ − . ξ ≈
23 (right) are plotted on the plane of slow variables. Phasetrajectories of the model plotted by dashed lines
For interpretation of the simulation results we shall utilize a scaled version of previously used slowvariables: x ∝ e cos ω, y ∝ − e sin ω. For exact relations between ( e, ω ) and ( x, y ), as well as the expression for small parameter ε , see AppendixA. Phase trajectories of several objects are plotted in Figure 22. Main sources of difference betweenanalytical phase portraits and numerical ones are non-zero eccentricity of Neptune ( e (cid:48) ≈ . . topological structure of phase portraits is reproduced in analytical model. Note, that in restricted three-body problem all solutions of Sessin and Ferraz-Mello (1984) on the same phase plane would be representedby concentric circles with e = const and ω changing linearly with time.The principal difference between our model and the one introduced by Sessin and Ferraz-Mello (1984)is the expansion of disturbing function past the first term in the Fourier series. Thus we can assert, that theregion of the complete phase space, where the second term influences the dynamics in a substantial way,is significantly large. Indeed about 10% of objects in our simulations deviated from circular trajectoriesin projection on the plane of slow variables. The rest correspond to very high or very low values of ξ , atwhich all trajectories in Q and Q , as well as the curve Γ separating them, in our model are likewisevery close to circles. Fig. 23
Resonant angle of 2011 UG vs time, demonstrating jumps between two different librating solutions in Q Some effects, that can be derived from the analytical model, are also observed in numerical simu-lations within the pool of selected KBOs. E.g. resonant angle of 2011 UG shows jumps between twodifferent librating regimes characteristic to motion in region Q (Fig. 23), while 2007 JJ demonstratesthe intermittent behaviour (Fig. 24) with the resonance angle constantly switching between libration andcirculation. Phase trajectory of 2007 JJ on the plane of slow variables is presented in Figure 25, showingthat changes in resonant angle behaviour are conditioned by secular trajectory crossing of the criticalcurve Γ . Similarly, analytical trajectory γ goes between regions of libration and circulation of resonantangle (Fig. 25). To compensate for angle ω precession on this kind of trajectories the modified resonantangle θ = ϕ + ω was used in previous plots . Same as ϕ , resonant angle θ circulates in Q and librates in Q and Q . Fig. 24
Resonant angle of 2007 JJ vs time. Intervals of circulation and libration colored green and gray respec-tively The resonant angle θ is the “standard” one for studies of first-order MMR (Murray and Dermott 2000). Ourchoice of the resonant angle ϕ was determined by the desire to write down the Hamiltonian of the system in theform (2), which seem to us most convenient for the analysis.2 Fig. 25
Phase trajectory of 2007 JJ alongside the analogous trajectory γ on the analytical phase portrait.Intervals of resonant angle circulation and libration on the left panel are colored green and gray respectively. Itis seen, that the green segments of the trajectory concentrate inside the region bound by Γ , while the gray onesmostly lie outside of it In this paper, using the averaging technique, we study Hamiltonian system that approximately describesthe dynamics of a three-body system in first-order MMR (within the restricted circular problem). Ourmodel incorporates harmonics of the Fourier series expansion of disturbing function up to the second order.Thus it can accurately describe the dynamics in that part of the phase space where first two harmonics havecomparable magnitudes, and where the well known integrable model of first-order MMR is not applicable.This is the region, from which chaos emerges.The nonintegrability of our model does not become an obstacle for a detailed analytical investigationof its properties. In particular, we have constructed bifurcation diagrams and phase portraits character-izing the long-term dynamics on different level sets defined by the system’s Hamiltonian and obtainedexpressions for probabilities of quasi-random transitions between different phase trajectories.The important question is the scope of the correctness of proposed model. We will try to answer it insubsequent studies, using Wisdom’s approach to investigate several first-order MMRs without truncatingthe averaged disturbing function. Nevertheless, even now we can note that the secular evolution of someKuiper belt objects qualitatively resembles what our simple model predicts.We hope also that our investigation outlines an approach which can be applied for analysis of similardegeneracy of the averaged disturbing function in the case of other MMRs (e.g., Sidorenko (2006)).
Appendix A: Constructing a model system, that reveals the origin of chaos in first-orderMMR
We shall confine ourselves to a case of exterior resonance p : ( p + 1) in restricted three-body problem.The interior resonance ( p + 1) : p can be reduced to the same model using similar approach. The distancebetween the two major bodies, i.e. a star and a planet, and the sum of their masses are taken here asunits of length and mass. The unit of time is chosen such that the orbital period of major bodies’ rotationabout barycenter is equal to 2 π . The mass of the planet µ is considered to be a small parameter of theproblem.Equations of motion for minor body (asteroid) in canonical form: d ( L, G, H ) dt = − ∂ K ∂ ( l, g, h ) , d ( l, g, h ) dt = ∂ K ∂ ( L, G, H ) , (28) where L, G, H, l, g, h are the Delaunay variables (Murray and Dermott 2000). They can be expressed interms of Keplerian elements a, e, i, Ω, ω as L = (cid:112) (1 − µ ) a, G = L (cid:112) − e , H = G cos i, g = ω, h = Ω. The last variable l is the asteroid’s mean anomaly.Hamiltonian K in (28) is K = − (1 − µ ) L − µR ( L, G, H, l, g, h − λ (cid:48) ) . (29)Here R is disturbing function in restricted circular three-body problem. Mean anomaly λ (cid:48) appearing in(29) depends linearly on time: λ (cid:48) = t + λ (cid:48) . Therefore it is convenient to use variable ˜ h = h − λ (cid:48) insteadof h , as it enables writing down the equations of motion in autonomous form as canonical equations withHamiltonian ˜ K = K ( L, G, H, l, g, ˜ h ) − H. We introduce resonant angle ¯ ϕ using the canonical transformation ( L, G, H, l, g, ˜ h ) → ( P ϕ , P g , P h , ¯ ϕ, ¯ g, ¯ h )defined by generating function S = ( p + 1) P ϕ l + (cid:2) P h + p (cid:0) P ϕ − P ∗ ϕ (cid:1)(cid:3) ˜ h + (cid:2) P g + ( p + 1) (cid:0) P ϕ − P ∗ ϕ (cid:1)(cid:3) g, where P ∗ ϕ = L ∗ / ( p + 1), while L ∗ = (cid:112) ( p + 1) /p is the value of L corresponding to exact p : ( p + 1) MMRin unperturbed problem ( µ = 0). The new variables are related to the old ones as follows: L = ∂S∂l = ( p + 1) P ϕ , G = ∂S∂g = P g + ( p + 1) (cid:0) P ϕ − P ∗ ϕ (cid:1) ,H = ∂S∂ ˜ h = P h + p (cid:0) P ϕ − P ∗ ϕ (cid:1) , ¯ ϕ = ∂S∂P ϕ = ( p + 1) l + p ˜ h + ( p + 1) g, ¯ g = ∂S∂P g = g, ¯ h = ∂S∂P h = ˜ h. The resonant angle can also be expressed in traditional form¯ ϕ = ( p + 1) λ − pλ (cid:48) − Ω. New Hamiltonian ˜ K = − (1 − µ ) p + 1) P ϕ − (cid:2) P h + p ( P ϕ − P ∗ ϕ ) (cid:3) − µR. The resonant case we are interested in corresponds to region R of phase space, selected by the condition (cid:12)(cid:12) pn (cid:48) − ( p + 1) n (cid:12)(cid:12) (cid:46) µ / . Here n (cid:48) = 1 and n are mean motions of planet and asteroid respectively. It is also true in R that (cid:12)(cid:12) P ϕ − P ∗ ϕ (cid:12)(cid:12) (cid:46) µ / . (30)In the resonant case, i.e. when the previous inequation holds, variables can be divided into fast , semi-fast and slow . Fast and semi-fast variables in R are ¯ h and ¯ ϕ respectively: d ¯ hdt ∼ , d ¯ ϕdt ∼ µ / . Slow variables, which vary with a rate of order µ , are P ϕ , P g , P h , and ¯ g .To study secular effects the averaging over the fast variable ¯ h is performed, which results in equationsof motion taking a canonical form with Hamiltonian¯ K ( P ϕ , P g , P h , ¯ ϕ, ¯ g ) = = 12 π ( p + 1) π ( p +1) (cid:90) ˜ K ( L ( P ϕ ) , G ( P g , P ϕ ) , H ( P h , P ϕ ) , l ( ¯ ϕ, ¯ g, ¯ h ) , ¯ g, ¯ h ) d ¯ h, where l ( ¯ ϕ, ¯ g, ¯ h ) = 1 k + 1 (cid:2) ¯ ϕ − ( k + 1)¯ g − k ¯ h (cid:3) . After such averaging the fast variable ¯ h vanishes, and the term “fast” is exempted. Thus in the rest ofthe paper we adopt name fast for denoting variables, which vary with the rate µ / , instead of referringto them as semi-fast, when the distinction of three different time scales was needed.Moreover, because there is no longer ¯ h in the Hamiltonian ¯ K , the conjugate momentum P h is a constantin considered approximation and can be treated as a parameter of the problem. Thus ¯ K is the Hamiltonianof a system with two degrees of freedom. Further instead of P h we shall use parameter σ = (cid:115) − P h L ∗ . Inequation e ≤ σ defines region S in phase space, to which the motion of the system is bound by theKozai-Lidov integral (Sidorenko et al. 2014).Next standard step in analysis of system’s dynamics in resonant region R is the scaling transformation(Arnold et al. 2006): ¯ τ = ¯ εt, ¯ Φ = ( P ∗ ϕ − P ϕ ) / ¯ ε, (31)where ¯ ε = µ / is a new small parameter. Using variables (31), the equations of motion can be rewrittenin a form of slow-fast system without loss of accuracy: d ¯ ϕd ¯ τ = χ ¯ Φ, d ¯ Φd ¯ τ = − ∂ ¯ W∂ ¯ ϕ ,dP g d ¯ τ = ¯ ε ∂ ¯ W∂ ¯ g , d ¯ gd ¯ τ = − ¯ ε ∂ ¯ W∂P g . (32)Here χ = 3 p / ( p + 1) / , ¯ W ( ¯ ϕ, ¯ g, P g ; σ ) = 12 π ( k + 1) π ( k +1) (cid:90) R (cid:16) L ∗ , P g , L ∗ (cid:112) − σ , l ( ¯ ϕ, ¯ g, ¯ h ) , ¯ g, ¯ h (cid:17) d ¯ h. For σ (cid:28) W ( ¯ ϕ, ¯ g, P g ; σ ) can be obtained as a series expansion (Murrayand Dermott 2000):¯ W ( ¯ ϕ, ¯ g, P g ; σ ) ≈ W + W e cos( ¯ ϕ − ¯ g )+ + e [ W + W cos 2( ¯ ϕ − ¯ g )] + i [ W + W cos 2 ¯ ϕ ] , (33)where e ≈ − P g L ∗ , i ≈ (cid:18) − P h L ∗ √ − e (cid:19) ≈ σ − e . (34)Expressions (34) are obtained taking into account resonance condition (30), and that the values e and i are limited by σ (cid:28)
1, leading to e, i (cid:28) W , W , W , W , W , W in (33) are calculated as follows: W = α b (0)1 / , W = α (cid:20) (2 p + 1) b ( p )1 / + α db ( p )1 / dα (cid:21) − δ p α ,W = α (cid:18) db (0)1 / dα + α d b (0)1 / dα (cid:19) ,W = α (cid:18)(cid:104) − p + 1) + 16( p + 1) (cid:105) b (2 p )1 / + α [8( p + 1) − db (2 p )1 / dα + α d b (2 p )1 / dα (cid:19) ,W = − W , W = α b (2 p +1)3 / . (35) Here α = ( p/ ( p + 1)) / , δ mn is the Kronecker delta, and b ( n )1 / ( α ), b ( n )3 / ( α ) are Laplace coefficients. Nu-merical values of coefficients (35) for several resonances are gathered in Table 1. Table 1
Numerical values of coefficients in series expansion of averaged disturbing function p = 1 p = 2 p = 3 p = 4 p = 5 W . . . . . W . . . . . W . . . . . W . . .
260 20 .
133 29 . W . . . . . Out of the whole region S we are interested in the part with small eccentricities. Thus we furtherassume e/σ (cid:46) σ , which leads to e (cid:46) σ . Taking into account (34) we can write down expression foraveraged disturbing function up to the terms of order σ :¯ W ( ¯ ϕ, ¯ g, P g ; σ ) ≈ W + W e cos( ¯ ϕ − ¯ g ) + σ [ W + W cos 2 ¯ ϕ ] . By introducing the variables¯ x = (cid:112) L ∗ − P g ) cos ¯ g, ¯ y = − (cid:112) L ∗ − P g ) sin ¯ g, the equations of motion (32) can be reduced to a Hamiltonian form d ¯ ϕd ¯ τ = ∂ ¯ Ξ∂ ¯ Φ , d ¯ Φd ¯ τ = − ∂ ¯ Ξ∂ ¯ ϕ ,d ¯ xd ¯ τ = ¯ ε ∂ ¯ Ξ∂ ¯ y , d ¯ yd ¯ τ = − ¯ ε ∂ ¯ Ξ∂ ¯ x with the Hamiltonian ¯ Ξ = χ ¯ Φ W √ L ∗ cos ¯ ϕ · ¯ x − W √ L ∗ sin ¯ ϕ · ¯ y + σ W cos 2 ¯ ϕ. The final rescaling of variables ϕ = − ¯ ϕ, Φ = − (cid:114) χσ W ¯ Φ,x = W σ W √ L ∗ ¯ x, y = W σ W √ L ∗ ¯ y,τ = σ (cid:112) χW ¯ τ , ε = W χ / L ∗ σ W / ¯ ε (36)results in the model system (1)–(2). Appendix B: Analytical expressions for integrals, emerging during the averaging over fastsubsystem’s period
Right-hand side parts of evolution equations (4), as well as adiabatic invariant formula (22), can beexpressed in terms of integrals I k,r = λ ∗ (cid:90) λ ∗ λ r dλ ( λ + 1) k (cid:112) ± ( λ − a )( λ − a )( λ − a )( λ − a ) , where k = 0 , , r = 0 ,
1, and integration limits λ ∗ , λ ∗ can be either real numbers or ±∞ . These integralscan be reduced to the linear combinations of elliptic integrals of the first, the second and the third kind: I , = c , K ( k ) ,I , = c , K ( k ) + c , Π ( h, k ) + ¯ c , Π (cid:0) ¯ h, k (cid:1) ,I , = g , K ( k ) + g , Π ( h, k ) + ¯ g , Π (cid:0) ¯ h, k (cid:1) ,I , = I , + c , K ( k ) + c , E ( k ) + c , Π ( h, k ) + ¯ c , Π (cid:0) ¯ h, k (cid:1) ,I , = g , K ( k ) + g , E ( k ) + g , Π ( h, k ) + ¯ g , Π (cid:0) ¯ h, k (cid:1) . (37)Further the formulae for coefficients c m,l , g m,l , moduli k and parameters h are gathered for all necessarycases.The case a j ∈ R ( j = 1 , a j are numbered in ascending order: a < a < a < a . Four different instances should be considered:A. Integration over ( a , a );B. Integration over ( a , a );C. Integration over ( a , a );D. Integration over ( −∞ , a ) (cid:83) ( a , + ∞ ). Instance A
Here λ ∗ = a and λ ∗ = a . We shall also use the following auxiliary parameters: A = 2 (cid:112) ( a − a )( a − a ) , C = 1 a − i , α = a − a a − a , α = ( a − a )( i − a )( a − a )( i − a ) . The value of elliptic integrals’ modulus in (37): k = (cid:114) ( a − a )( a − a )( a − a )( a − a ) , and parameter h = α .In order to find coefficients in (37), we considered such linear combinations of I k,r , that are reducedto integral (253.39) from (Byrd and Friedman 1954): I , + iI , = a (cid:82) a dλ ( λ − i ) √ − ( λ − a )( λ − a )( λ − a )( λ − a ) ,I , − iI , − I , = − a (cid:82) a dλ ( λ − i ) √ − ( λ − a )( λ − a )( λ − a )( λ − a ) . (38)After some simple calculations, one can find: g , = Re C , c , = Im C , C = A C α α ,g , = A C (cid:16) − α α (cid:17) , c , = − ig , ,c , = − Re C , g , = Im C , C = A C α (cid:20) α + ( α − α ) α − (cid:21) ,c , = − Re C , g , = Im C , C = A C ( α − α ) α ( α − k − α ) ,c , = − A C ( α − α )4 α (cid:104) α + (2 α k +2 α − α − k )( α − α )2( α − k − α ) (cid:105) ,g , = ic , . (39)Using (253.00) from Byrd and Friedman (1954) we also obtain: c , = A . (40) Instance B
For integrals over the interval ( a , a ) the only difference is the parameters α , α and modulus k . Thusin (39) the values α = a − a a − a , α = ( a − a )( i − a )( a − a )( i − a ) , k = (cid:114) ( a − a )( a − a )( a − a )( a − a )should be used. Instance C
For integrals over the interval ( a , a ) in (39): C = 1 a − i , α = a − a a − a , α = ( a − a )( i − a )( a − a )( i − a ) . The modulus k and parameter A are the same as for interval ( a , a ). Note:
The equality of k and c , in instances A and C means, that when the potential (2) has twominima, the periods of librations T about them on the same energy level are equal. This also holds truefor local minima corresponding to instances B and D , as the substitution λ = tan[( ϕ − ˜ ϕ ) /
2] always allowsto get rid of semi-infinite intervals of integration. Moreover, the averaged values of cos in (18) are also thesame for librations around two local minima. The averaged values of sin in (18) for librations about twolocal minima differ by π , and values of adiabatic invariant (22) differ by y/ √ Instance D
For integrals over two semi-infinite intervals the values of parameters in (39) are C = 1 a − i , α = a − a a − a , α = ( a − a )( i − a )( a − a )( i − a ) . The modulus k is the same as for integration over ( a , a ), and A is the same as for ( a , a ).The case a , a ∈ R ( a < a ), a , a ∈ C ( a = ¯ a )Here there are only two instances:A. Integration over ( a , a );B. Integration over ( −∞ , a ) (cid:83) ( a , + ∞ ). Instance A
Here the following auxiliary parameters in expressions for c m,l and g m,l will be used: α = ( a A − a A ) − ( A − A ) i ( a A + a A ) − ( A + A ) i , α = A − A A + A ,C = A + A ( A a − A a ) − ( A − A ) i , where A = (cid:112) ( a − a ) + b , A = (cid:112) ( a − a ) + b , A = 1 / √ A A ,a = Re a = Re a , b = Im a = − Im a . Using (38) and formulae (259.04), (341.01)–(341.04) from (Byrd and Friedman 1954), we obtain c , = 2 A , g , = Re C , c , = Im C , C = 2 A C α ,c , = − i ( α − α )(1 − α ) A C , g , = ic , ,c , = − Re C , g , = Im C , C = A C (cid:104) α + ( α − α ) α − (cid:105) ,c , = − Re C , g , = Im C , C = A C α ( α − α ) (1 − α )( k + α k (cid:48) ) ,c , = A C α − α )( α − (cid:20) α + ( α − α ) [ α (2 k − − k ] ( α − k + α k (cid:48) ) (cid:21) , g , = ic , . (41)Modulus and parameter of elliptic integrals are calculated as follows: k = (cid:115) ( a − a ) − ( A − A ) A A , h = α α − . Instance B
Expressions (41) stay the same. Auxiliary parameters are α = ( a A + a A ) − ( A + A ) i ( a A − a A ) + ( A − A ) i , α = A + A A − A ,C = A − A ( A a + A a ) − ( A + A ) i . The modulus k of elliptic integrals is also different: k = (cid:115) ( A + A ) − ( a − a ) A A . The case a j ∈ C ( j = 1 , a j as p , , and imaginary parts as b , : a , = p ± b , a , = p ± b . Without loss of generality let us assume, that p < p , b >
0, and b > p = p , then also b < b ).In expressions for c m,l and g m,l the following auxiliary parameters will be used A = (cid:113) ( p − p ) + ( b − b ) , A = (cid:113) ( p − p ) + ( b + b ) , A = 2 / ( A + A ) ,g = (cid:115) b − ( A − A ) ( A + A ) − b , α = b + g ( p − i ) p − b g − i , C = 1 b + g ( p − i ) . Similar to previous case we find: c , = 2 A , c , = Im C , g , = Re C , C = α (1+ g α )1+ α A C ,g , = α ( α − g ) A C α , c , = − ig , , c , = − Re C , g , = Im C ,C = A C (cid:104) g + g ( α − g )1+ α + ( α − g ) (1+ α )( α + k (cid:48) ) (cid:16) k (cid:48) +2 α − α k α − k (cid:48) (cid:17)(cid:105) ,c , = Re C , g , = − Im C , C = A C ( α − g ) α ( α +1)( α + k (cid:48) ) ,c , = − A C α − g ) α ( α +1) (cid:104) g + ( α − g )(2 k (cid:48) +2 α − α k )( α +1)( α + k (cid:48) ) (cid:105) , g , = ic , . Modulus and parameter of elliptic integrals: k = 2 √ A A A + A , h = α + 1 . It should be noted, that (259.04) in (Byrd and Friedman 1954) contain a typo – the multiplier g (in authors’notation) is missing.9 Acknowledgements
The work was supported by the Presidium of the Russian Academy of Sciences (Program 28“Space: investigations of the fundamental processes and their interrelationships”). We are grateful to A.I.Neishtadt,A.Correia, A.Morbidelli and J.Wisdom for useful discussions. We would also like to thank D.A.Pritykin for proof-reading the manuscript.