Secondary resonances and the boundary of effective stability of Trojan motions
PP´aez & Efthymiopoulos (2017) Submitted to CMDA
Secondary resonances and the boundary of effective stabilityof Trojan motions
Roc´ıo Isabel P´aez and Christos Efthymiopoulos Research Center for Astronomy and Applied Mathematics, Academy of Athens
Abstract:
One of the most interesting features in the libration domain of co-orbital motions is the existenceof secondary resonances. For some combinations of physical parameters, these resonances occupya large fraction of the domain of stability and rule the dynamics within the stable tadpole region.In this work, we present an application of a recently introduced ‘basic Hamiltonian model’ H b forTrojan dynamics [33], [35]: we show that the inner border of the secondary resonance of lowermostorder, as defined by H b , provides a good estimation of the region in phase-space for which the orbitsremain regular regardless the orbital parameters of the system. The computation of this boundary isstraightforward by combining a resonant normal form calculation in conjunction with an ‘asymmetricexpansion’ of the Hamiltonian around the libration points, which speeds up convergence. Applicationsto the determination of the effective stability domain for exoplanetary Trojans (planet-sized objectsor asteroids) which may accompany giant exoplanets are discussed. Despite the theoretical possibility of the existence of Trojan exoplanets ([18], [1], [3]), no such body hasbeen identified so far in exoplanet surveys. This lack of identification may reflect formation constrains,constrains to detectability ([17], [4], [19], [20], [21]), or it may simply be due to stability reasons. Inthis framework, the question of ‘effective stability’, i.e. stability of the orbit of a Trojan body fortimes as long as a considerable fraction of the age of the hosting system, comes to the surface. Thequestion of effective stability has been addressed nearly exhaustively in the case of Trojan asteroidsin our Solar System (see, for example, [27], [22], [15], [39], [37], [23], [5], [24]) from both numericaland analytical approaches, but only scarcely in the case of exoplanetary systems (see [30], [12], [38],[9]). One main reason for the scarcity of results in this latter case is the vast volume of parameterspace to be investigated, in conjunction with the multi-body nature of the problem: to determinethe long-term stability of Trojan motions becomes essentially a problem of secular dynamics with asmany degrees of freedom as the number of planets in the system under consideration. Any attemptto face the problem other than numerical simulation clearly requires a simplification of the dynamicalmodel, without this leading to oversimplified conclusions regarding the long-term orbital stability.In the present work, we discuss a key property of the dynamics induced by secondary resonances in the domain of Trojan motions, which in addition to its own proper interest, can serve also thepurpose of obtaining a simple analytical estimate of the effective stability boundary of Trojan motionsin hypothetical exoplanetary systems. Our analysis of the resonant dynamics stems from a set ofconsiderations or assumptions, whose validity can be most easily judged by comparison with someresults and figures of a previous work of ours ([33]) as follows:1) In [33] we provided a formalism of the problem of the dynamics of Trojan bodies in the Hamiltoniancontext, which recovers all essential features as discovered in previous literature ([10], [11], [28], [29]);preceding works, however, focus mostly on a direct investigation of the equations of motion, averagedor not with respect to short period terms. In our works we stressed, instead, a main advantage of thenew formalism, namely the allowance to recruit the full machinery of Hamiltonian methods in orderto better analyze the problem under study.2) We investigated the features of Trojan dynamics which hold under three physically relevantassumptions: i) that the motions of all planets, including the Trojan body, are close to planar, ii)that the Trojan body is small enough to be considered as test particle, in agreement with formationscenaria which suggest that exo-Trojans should be at most Mars-sized objects ([1]), and iii) that the a r X i v : . [ a s t r o - ph . E P ] D ec ´aez & Efthymiopoulos (2017) Submitted to CMDA secular dynamics of the hosting system is such that the eccentricity vector of the primary companionof the Trojan body undergoes circulation with a nearly constant frequency g (cid:48) , and has a length whichundergoes variations around some non-zero value e .Under assumptions (i) to (iii), we find that the Hamiltonian of motion of the Trojan body, averagedover short period terms for the motions of the remaining planets, can be decomposed in the form H = H b + H sec where: H b , called the ‘basic model’, describes short period and synodic motions, andyields a constant proper eccentricity for the Trojan body, and H sec contains all remaining secularperturbations. Furthermore H b has a universal form, i.e., it suffices to redefine the physical meaningof the angular canonical variables, to keep its form unaltered in the whole hierarchy of restrictedproblems (circular, elliptic, secular with more than one perturbing planets).3) The decomposition H = H b + H sec leads to a specific physical understanding of the dynamicswhen the primary planet has a mass in the giant planet range. In this case, the three timescalesrelated to the short-period, synodic and secular motions have a separation by about one or less orderof magnitude from each other. Then, due to the specific features of H b described above, we arriveat the following key remark: the model H b produces, in phase space, a set of secondary resonancescorresponding to commensurabilities between the frequencies of the short-period and synodic mo-tions ([13]). It is easy to see that these are the only secondary resonances which occupy a non-zerovolume in phase space in the whole hierarchy of restricted problems that one could use as dynamicalmodels for the Trojan body. However, there exists a modulation effect ([2]) due to the influence of H sec on this set of secondary resonances: the separatrices pulsate slowly (with one or more secular frequen-cies) and, as a result, in the ‘domain of uncertainty’ ([31]) created by such pulsations, the motionsbecome chaotic. Such an effect is possible to visualize already in the Elliptic Restricted Three-BodyProblem, namely the simplest model with non trivial H sec . The reader is refered to Figures 5 to 15of [33] which show in detail the statements below, by exemplifying the outcome of the modulationeffects for the secondary resonances 1:5 up to 1:12, when the modulus of the eccentricity vector e of the primary companion varies from e = 0 to just a moderate value e = 0 .
1. By inspecting thestability maps in the space of the Trojan body’s proper elements, one sees that, for e slightly largerthan zero, the separatrix pulsation for the secondary resonances becomes large enough so as to wipeout nearly completely the domain of stable motions occupied by such resonances. As a result, theonly remaining stable motions are those in a inner (closer to the libration center) domain devoid ofsecondary resonances. In fact, as found in many works (e.g. [37], [26]) there can still be resonancesinvolving one or more secular frequencies which penetrate this innermost stability domain. However,since these resonances are thin and typically do not overlap, they can only induce a very slow chaoticdiffusion of the Arnold type, which can be neglected for all practical purposes. Hence, the innermostdomain, devoid of the secondary resonances of H b , meets all criteria of effective stability, and, indeed,stability maps indicate the robustness of this domain against variations of the orbital parameters ofthe Trojan body. Stemming from remarks (1) to (3) above, we propose below a practical method to define the effectivestability domain of Trojan motions. This is based on the following steps:
Step 1: analyze a given system where hypothetical Trojan bodies are sought for and compute theHamiltonian H b , Step 2: identify the largest in size (typically lowest in order) secondary resonance of H b for givenparameter values, Step 3: compute a resonant normal form and evaluate the theoretical separatrices of the identifiedsecondary resonance,
Step 4: assume that all stable domains of resonant motions extending beyond the innermost (closestto the libration center) branch of the theoretical separatrices were wiped out by secular modulationeffects.Then, the locus S formed by the intersection of the family of all computed innermost theoreticalseparatrices, along the dominant secondary resonance, with any chosen (with respect to phases) plane P P of Trojan proper elements, yields the boundary of the effectively stability domain in the plane
P P . The above computation is fast and straightforward to perform with modern computer algebraprograms, and thus competitive to large grid computations of effective stability maps. In our ownimplementation we use the normal form method adopted in [34], [35].In the rest of the paper, we discuss both the dynamical role of the secondary resonances in delimit-ing the main domain of effective stability as well as our particular analytical method of computing theborder of this domain. The structure of the paper is as follows: In Section 2, we review the derivationand features of the Hamiltonian H b , as well as our way to expand H b in a form form suitable forresonant normal form computations. A novel feature is the adoption of an ‘asymmetric expansion’which improves convergence. Section 3 explains in detail the realization of the Steps 1-4 , in partic-ular the computation of the theoretical separatrices of the dominant secondary resonance and theirsuperposition to stability maps in the space of proper elements. Section 4 contains the main results:(a) we provide numerical evidence, based on stability maps, of how the separatrices of the secondaryresonances of the ’basic model’ H b act as delimiters of the effective stability domain; (b) we use ananalytical method to estimate this boundary; (c) we discuss the robustness of the present approachagainst changing the model’s parameters (masses and eccentricities), as well as when considering, inthe numerical integrations, the full three-body problem instead of the ERTBP. Section 5 summarizesour main conclusions. H b and its asymmetric expansion H b In [33], a Hamiltonian formulation was provided for the Trojan motion which applies to the planarElliptic Restricted Three-Body Problem (ERTBP) with a central mass, a primary perturber or simply‘primary’, and the Trojan test particle, or when S additional perturbing bodies are present but farfrom MMRs, the so-called ’Restricted Multi-Planet Problem’ (RMPP)). The Hamiltonian reads H = H b ( Y f , φ f , u, v, Y p ; µ, e (cid:48) ) + H sec ( Y f , φ f , u, v, Y p , φ, Y , φ , . . . , Y S , φ S ) . (1)In Eq. (1), the variables ( φ f , Y f ), ( u, v ) and ( φ, Y p ) are pairs of action-angle variables, whose definitionstems from Delaunay-like variables following a sequence of four consecutive canonical transformations(see Appendix A). In particular, ( Y f , φ f ) are action-angle variables describing the fast degree offreedom, of frequency ω f ≡ ˙ φ f = 1 − µ + g (cid:48) + . . . , (2)where g (cid:48) is fundamental frequency of precession of the primary’s perihelion. The pair ( u, v ) describethe particle’s synodic librations, u (cid:39) λ − λ (cid:48) − π/ v (cid:39) √ a −
1, with λ , λ (cid:48) the mean longitudes ofthe test particle and of the primary, a the particle’s major semi-axis, and a (cid:48) = 1. The associatedfrequency at the libration center is ω s ≡ ˙ φ s = − (cid:114) µ . . . . (3)Finally, the secular motion of the test particle’s eccentricity vector ( e cos( ω − ω (cid:48) ) , e sin( ω − ω (cid:48) )), where e is the eccentricity and ω , ω (cid:48) are the arguments of the perihelion of the test particle and the primaryrespectively, is described by a circulation around the forced equilibrium point, given in our variablesby a set of action angle variables ( Y p , φ ). The associated secular frequency is g ≡ ˙ φ = 278 µ − g (cid:48) + . . . . (4)We call the term H b in the Hamiltonian of Eq. (1) the ‘basic Hamiltonian model’ for Trojan motionsin the 1:1 MMR. Its detailed form is given in the Supplementary Online Material of [33]. We find H b = − v ) − v + (1 + g (cid:48) ) Y f − g (cid:48) Y p − µ F (0) ( u, φ f , v, Y f − Y p ; e (cid:48) ) . (5) The physical parameters entering into H b are i) the mass parameter µ = m (cid:48) m (cid:48) + M , where M is the massof the central mass and m (cid:48) the mass of the primary, ii) the mean value of the length of the eccentricityvector of the heliocentric orbit of the primary perturber, e (cid:48) . In the ERTBP, one has simply e (cid:48) = e (cid:48) , g (cid:48) = 0 (implying also ω (cid:48) ≡ const). However, the form of H b remains the same in both the ERTBPand the RMPP. In particular, the angle φ is defined via a ‘shift transformation’ depending only on therelative difference ∆ ω = ω − ω (cid:48) . Physically, the secular dynamics induced under H b appears the samein the ERTBP and in the RMPP, when, in the latter case, it is viewed in apsidal co-rotation withthe primary. Furthermore, since the angle φ is ignorable in H b , the action variable Y p is an integralof the basic Hamiltonian. Then, the ERTBP and the RMPP are diversified only by their differentform of the functions H sec . In particular, in the RMPP case H sec contains also pairs of action anglevariables associated with the secular precessions of the S additional bodies, while in the ERTBP itcontains only the angle φ associated with the secular precession of the test particle. Finally, H sec disappears all together in the circular RTBP. Thus, H b becomes the exact Hamiltonian in this case.Note, however, that in the ERTBP H b is not equal to the ensemble of all terms independent of e (cid:48) .The basic model H b represents a drastic reduction of the number of degrees of freedom withrespect to the original problem. In the sequel, we will focus on one particular feature of H b , namelythe presence of secondary resonances , which correspond to commensurability relations between ω f and ω s . In particular, we focus on the role of these resonances in practically determining the boundary ofthe effective domain of stability for the Trojan motions. The resonant normal form computed in Section 3 below provides a model for studying the dynamicswithin or near a secondary resonance of the form: m f ω f + m s ω s = 0 . (6)A non-resonant normal form for the model H b , allows to find the location of secondary resonances ina space of suitably defined proper elements for the Trojan body (see [35]). However, the non-resonantnormal form does not allow to compute the local phase portrait, i.e., the separatrices associated witheach resonance. Furthermore, all series expansions which are polynomial in the variables u, v exhibitpoor convergence, a fact associated with the singularity (collision with the primary) at u = − π/
3. Inorder to deal with this problem, a partially expanded version of the H b can be used [34], in which allthe powers of the quantity β ( τ ) = √ − τ (with τ = u + π/
3) are kept unexpanded. This leads toa Hamiltonian of the form H b ( v, Y , τ, φ f , Y p ) = − v + ∞ (cid:88) i =0 ( − i − ( i + 1) v i Y + Y p + µ (cid:88) m ,m ,m k ,k ,k ,j a m ,m ,m ,k ,k ,j e (cid:48) k v m cos k ( τ ) sin k ( τ ) Y m cos m φ f sin m φ f β j ( τ ) , (7)where τ = u + π/ Y = Y f − Y p , and a m ,m ,m ,k ,k ,j are rational numbers.The librations in τ (or u ) are represented in terms of the synodic angle variable φ s , i.e., the phaseof the synodic libration. The computation of a resonant normal form requires to explicitly Fourierexpand the terms of H b in both angles φ f and φ s . Although the Hamiltonian (7) represents a Fourierexpansion for the fast d.o.f. (angle φ f ), there still remain the powers of β that must be expandedin powers of u in order to obtain a complete Fourier expansion in the angle φ s as well. Due to thesingularity at τ = 0 (or u = − π/ β ( τ ) N = − τ ) N/ ,with N ∈ N , around a certain τ is convergent only in the domain D τ ,δ centered at τ and of radius δ = Min { τ , π − τ } . The most common approach consists of Taylor expansions around the librationequilibrium point, located at τ = π , for L4, or τ = π for L5. The corresponding δ in this case is π . One finds that many Trojan orbits, and important secondary resonances, may cross this domain.In such cases, the resonant normal form construction is obstructed by the poor convergence of theoriginal Hamiltonian expansion. In order to face this problem, we find a different polynomial representation of the Hamiltonian H b in the variables ( u, v ) by performing an asymmetric expansion, i.e. expansion around a non-equilibrium point τ (cid:54) = π/
3, selected to be further away from the singularity but close enough to thelibration point, so that a re-ordering of the expansion in powers of u yields a negligible term linear in u (since u = 0 represents the equilibrium point of H b ). Here we choose τ = π . In this case, we obtaina polynomial expansion of the Hamiltonian in powers of the quantity ( τ − π/ u . It is immediate to see that any finite truncationof this expression yields a different polynomial than the one obtained by a finite truncation of thedirect Taylor expansion around τ = π/
3. However, the new expression better represents the quantities β ( τ ) N in a domain extended up to τ ∼ π . We call the expansion around τ = π asymmetric , whilethe one around τ = π symmetric .Figure 1 shows the benefits of the asymmetric expansion when compared to the symmetric one.We consider the functions B ( τ ) = cos τβ ( τ ) = cos τ (2 − τ ) / , B ( τ ) = cos τβ ( τ ) = cos τ (2 − τ ) / ,B ( τ ) = cos τβ ( τ ) = cos τ (2 − τ ) / , B ( τ ) = cos τβ ( τ ) = cos τ (2 − τ ) / , (8)which represent the most common terms in powers of β ( τ ) appearing in Eq. (7). The symmetricTaylor expansion of B , B , B and B around τ = π/ B M,π/ ( u ) = B M ( π/
3) + B (1) M ( π/ u + 12 B (2) M ( π/ u + 16 B (3) M ( π/ u + . . . , (9)where B ( n ) M ( π/
3) is the n -th derivative of the function B M , evaluated at u = π/ M = 1 , , , u = τ − π/
3. On the other hand, the asymmetric Taylor expansions of the same functions around τ = π/ B M,π/ ( u ) = B M ( π/
2) + B (1) M ( π/
2) ( u − π B (2) M ( π/
2) ( u − π + 16 B (3) M ( π/
2) ( u − π + . . . . (10)Fig. 1 compares the graphs of the original functions B M ( u ) (pink) with the two correspondingexpansions B M,π/ ( u ) (symmetric, blue) and B M,π/ ( u ) (asymmetric, green) up to order 10 in u . Wesee that both expansions provide a good representation of the original function up to a certain extentin u , but for increasing values of u , the asymmetric expansions are more accurate than the symmetricones in a domain extending to higher values of u . The improvement in accuracy is more notoriousas M increases. In fact, we find that the polynomial approximations to H b found by the asymmetricexpansion is accurate up to u ∼ u ∼
0, the symmetric expansion is marginallymore accurate than the asymmetric one. This yields a slight shift of the equilibrium point of H b withrespect to u = 0, typically of about ∼ − rad, i.e. practically negligible.Similar results are found for the asymmetric and symmetric expansions of the functions sin τβ ( τ ) ap-pearing in H b . Finally, the asymmetric expansions for both types of functions can be easily performedby a closed set of formulas, given in the Appendix B. The construction of the asymmetrically expanded H b consists of two steps: i) replacement of theexpansions for the variables ( u, v ) in Eq. (7), and ii) transformation to action-angle variables.In order to replace the polynomial truncations for the functions of Eq. (8) into Eq. (7), we adoptthe asymmetric expansion (10), using the formulæ provided in the Appendix B. Regarding v , itis enough to consider the Taylor expansion of H b with respect to v around zero. The maximum B M ( u ) (pink), with the expansions B M,π/ ( u ) (blue)and B M,π/ ( u ) (green), up to order 10 in u , for M = 1 , , , u ∈ [0 . , truncation order is determined in terms of a ‘book-keeping parameter’ ([8]; see below). After thesereplacements, the Hamiltonian H b takes the form H b ( v, Y , u, φ f , Y p ) = Y p + (cid:88) m ,m ,m ,m a m ,m ,m ,m v m u m ( √Y ) m cos( m φ f ) , (11)where the real coefficients a m ,m ,m ,m depend on the parameters µ and e (cid:48) (or simply e (cid:48) in theERTBP).Next, we diagonalize H b in order to obtain a harmonic oscillator quadratic part for the synodicdegree of freedom. Diagonalization is performed by the linear canonical transformation ( u, v ) → ( U, V )defined by the set of formulas: (cid:18) uv (cid:19) = 1 (cid:112) Det( E ) ( E · B ) (cid:18) UV (cid:19) , B = (cid:32) √ − i √ − i √ √ (cid:33) , (12)where E is a 2 × e , associated with the eigenvalues λ , = ± ω s of the matrix M M = (cid:18) a (1 , , , a (2 , , , − a (0 , , , − a (1 , , , (cid:19) . (13)From the variables U and V we then pass to the action-angle variables ( Y , φ f ) with U = (cid:112) Y s sin φ s , V = (cid:112) Y s cos φ s . (14)The last step corresponds to a re-organization of the terms of the Hamiltonian, according to a book-keeping parameter [8]. This is a parameter with numerical value equal to (cid:15) = 1. To every termin the Hamiltonian (11), we asign a power of (cid:15) indicating the order of the normalization at whichthe term will be treated. Thus, coefficients with powers of (cid:15) propagate throughout the series at allnormalization steps, helping to organize the terms in different orders of smallness. Regarding theoriginal Hamiltonian, we adopt the following book-keeping rule: Rule 3.1
To every monomial of the type c ( k ,k ,k ,k ) ( (cid:112) Y s ) k ( √Y ) k cossin ( k φ s + k φ f ) , assign a book-keeping coefficient (cid:15) r ( k ,k ,k ) , where the exponent r ( k , k , k ) is given by r ( k , k , k ) = (cid:40) Max(0 , k + k − if k = 0Max(0 , k + k −
2) + 1 if k (cid:54) = 0 . This book-keeping rule ensures also that the terms of zero-th order in (cid:15) are linear in Y and Y s . TheHamiltonian now takes the form: H b ( Y s , Y , φ s , φ f , Y p ) = Y p + ω s Y s + ω f Y + r max (cid:88) r =1 c ( k ,k ,k ,k ) (cid:15) r ( (cid:112) Y s ) k ( √Y ) k cossin ( k φ s + k φ f ) . (15)From the canonical transformation in Eq. (14), it is straighforward to check that the harmonicsof the angles φ f and φ s have the same parity as the powers of the corresponding functions in thevariables √Y and √ Y s . The resonant normalization of the Hamiltonian (15) consists of a sequence of near-identity canonicaltransformations, in ascending powers of the book-keeping parameter (cid:15) , aiming to eliminate from theHamiltonian the trigonometric dependence on the angles in any linear combination other than the onewhich corresponds to the selected secondary resonance (Eq. 6). The resulting normal form includes,besides terms depending just on the actions, also terms of the form b ( p ( r ) ) e i( k · q ( r ) ) . (16)Here, q (0) = ( φ f , φ s ), p (0) = ( Y , Y s ), and the superscript ( r ) indicates the variables found after r consecutive near-identity normalizing transformations of ( q (0) , p (0) ). Also, k = ( k , k ) belongs to theset M called the resonant module (Eq. 17 below). The terms (16) allow to determine the theoreticalseparatrices of the secondary resonance via the process described in subsection 3.3 below.The general recursive resonant normalization algorithm is defined as follows: Let m , m be twointegers marking the secondary resonance m m ≈ ω f ω s . The resonant module M is the set of integervectors defined by M = { k = ( k , k ) : k m + k m = 0 } , (17)where (cid:80) i =1 | m i | (cid:54) = 0.Let us assume that the Hamiltonian is in normal form up to order r in the book-keeping parameter,i.e. H = Z + (cid:15) Z + . . . + (cid:15) r Z r + (cid:15) r +1 H ( r ) r +1 + (cid:15) r +2 H ( r ) r +2 + . . . . (18)From the terms of order (cid:15) r +1 , in the Fourier expansion, H ( r ) r +1 = (cid:88) k b ( p ( r ) ) e i( k · q ( r ) ) , (19)where we isolate the terms that we want to eliminate in the present step, denoted by ∗ H ( r ) r +1 = (cid:88) k / ∈M b ( p ( r ) ) e i( k · q ( r ) ) . (20)The homological equation (cid:15) r +1 ∗ H ( r ) r +1 + {Z , χ r +1 } = 0 (21) has the solution χ r +1 = (cid:15) r +1 (cid:88) k / ∈M b ( p ( r ) )i ( k · ω ) e i( k · q ( r ) ) , (22)with ω = ( ω f , ω s ).Having the expression of the generating function, we compute the transformed Hamiltonian H ( r +1) = exp( L χ r +1 ) H ( r ) , (23)where exp (cid:16) L χ (cid:17) · = I · +( L χ · ) + 12 ( L χ · ) + . . . . (24)and the Lie operator L χ ≡ {· , χ } ( {· , ·} denotes the Poisson bracket).By construction, the Hamiltonian in Eq. (23) is in normal form up to order (cid:15) r +1 , i.e. H = Z + (cid:15) Z + . . . + (cid:15) r Z r + (cid:15) r +1 Z r +1 + (cid:15) r +2 H ( r ) r +2 + (cid:15) r +3 H ( r ) r +3 + . . . . (25) Let us consider the function H b given in Eq. (15) as the starting Hamiltonian H (0) b of the normalizingscheme. We apply the normalizing scheme presented above, up to a maximum normalization order R in (cid:15) . In the examples that follow, the maximum normalization order examined was R = 22. However,since the resonant normal form series are asymptotic, depending on the parameters and resonanceconsidered, the optimal normalization order (yielding the minimum remainder as computed e.g. in [7])varies, yielding optimal orders between R = 14 and R = 20.Let H ( R ) b be the final normalized Hamiltonian. According to Eq. (16), the form of H ( R ) b is givenby H ( R ) b = R (cid:88) r =0( k f ,k s ) ∈M (cid:15) r b ( Y ( R ) , Y ( R ) s ) e i( k f φ ( R ) f + k s φ ( R ) s ) . (26)If we replace the book-keeping parameter (cid:15) for its value equal to 1, we recover the final normal form,depending on the actions and the angles through the combination, H ( R ) b = (cid:88) ( k f ,k s ) ∈M c ( d f ,d s ,k f ,k s ) (cid:112) Y ( R ) d f (cid:113) Y ( R ) s d s e i( k f φ ( R ) f + k s φ ( R ) s ) , (27)where the pairs ( d f , k f ) and ( d s , k s ) have the same parity, and the values of the Fourier wavenumbersare bounded by | k f | ≤ d f and | k s | ≤ d s . The integers ( d f , d s ) are limited by the value of R , throughthe book-keeping Rule 3.1.We define the quantity Ψ = m Y ( R ) + m Y ( R ) s (28)as a resonant integral of the normal form H ( R ) b , where m and m are the integers that define theresonant module M in Eq. (17). Considering Eq. (16), it is straightforward to prove that L H ( R ) b Ψ = { H ( R ) b , Ψ } = 0 , (29)i.e. Ψ is a formal integral of H ( R ) b .By considering the transformation C ( R ) , C ( R ) = ϕ (1) ◦ ϕ (2) ◦ . . . ◦ ϕ ( R − ◦ ϕ ( R ) , (30)where ϕ ( r ) = exp ( L χ r ) ( Y ( r ) , Y ( r ) s , φ ( r ) f , φ ( r ) s ) (31) we can represent the resonant integral in terms of the original variables ( Y (0) , Y (0) s , φ (0) f , φ (0) s ), viaΨ( Y (0) , Y (0) s , φ (0) f , φ (0) s ) = Ψ (cid:16) C ( R ) ( Y ( R ) , Y ( R ) s , φ ( R ) f , φ ( R ) s ) (cid:17) . (32)Finally, applying the inverse transformations to those of Eqs. (14) and (12), we are able to expressthe resonant integral in (32) as function of the variables used in Eq. (11)Ψ ≡ Ψ( v, Y , u, φ f ) . (33)Having arrived at a final expression for the resonant integral Ψ in terms of the original canonicalvariables, we can compute the form of the theoretical separatrices of the corresponding secondaryresonance in any suitably defined surface of section of the Hamiltonian H b . In the numerical resultsbelow, we adopt a section of the form φ f = φ f , as well as a constant value of the energy E = H b ,the equation E = H b ( v, Y , u, φ f ) can be solved for Y . Substitution to (33) yields then the resonantintegral on the surface of section as a function of u and v only, viz.Ψ ≡ Ψ (cid:16) v, Y ( u, v ; E, φ f ) , u, φ f (cid:17) . (34)The theoretical phase portrait is now obtained by the level curves of Eq. (34). Figure 2, left panel,summarizes the main features of the theoretical phase portrait. In particular, the stable periodic orbitof the secondary resonance is represented by the points of extremum of the level set of Ψ, while theunstable periodic orbit corresponds to the minimax (saddle) points of the level set of Ψ. The levelcurves with Ψ = Ψ nmx , where Ψ nmx is the value of the resonant integral at the saddle points, are thecurves representing the theoretical separatrices of the secondary resonance. We present below numerical results based on the computation of stability maps for selected valuesof the parameters µ and e (cid:48) , characterized by the presence of conspicuous secondary resonances ofthe Hamiltonian H b . The stability maps are given in color scale of the values of the Fast LyapunovIndicator ([14]), for orbits with initial conditions labeled in terms of two quantities (∆ u, e p ). Thesequantities also serve as proper elements, i.e. quasi integrals of motion, for the subset of all regularorbits in every stability map. Working on fixed surfaces of section φ f = − π/
3, the relation betweeninitial conditions ( u, v, Y ) and (∆ u, e p ) is given by the relations Y = e p, , ∆ u = u − u , where u is the point of intersection of the short-period orbit around L4 with the surface of section (see [33]for analytical expressions), and v = B ∆ u , for fixed parameters B (depending on µ ) selected in sucha way that the straight line v = B ( u − u ) in the surface of section passes right through one of theislands of the secondary resonance chain. The half-witdh of the libration in u as a function of ∆ u , B , µ , e p and e (cid:48) , reads D p = (cid:34) B / µ (cid:0) / e (cid:48) /
16 + 129 e p / (cid:1) µ (cid:0) / e (cid:48) /
16 + 129 e p / (cid:1) (cid:35) / ∆ u + O (∆ u ) . (35)The values of B used in the various stability maps below are given explicitly in the caption of eachfigure.We can now superpose the theoretical computation of the phase portrait of the secondary reson-ance to the numerical results found in the stability maps. For given parameters µ, e (cid:48) , B , and choosingone value of the energy E , one obtains the resonant integral (34) as a function of u only. An exampleis shown in the right panel of Fig. 2. The value u = u marks the position of local maximum of theresonant integral Ψ along the line v = B ( u − u ). This corresponds a central locus passing approx-imately through the middle of the resonant domain along the corresponding secondary resonance. u, v ) for a surface of sectionof the H b . The central blue dot represents the location of a stable periodic orbit, whoseco-ordinate in equal to u = u res . At this point, the resonant integral Ψ presents a globalextremum. Additional quasi-periodic orbits inside the island of stability are labeled with thecorresponding values of Ψ, i.e. Ψ ∗ , Ψ ∗ , Ψ ∗ , Ψ mnx , accomplishing Ψ ∗ > Ψ ∗ > Ψ ∗ > Ψ mnx . The value Ψ mnx represents a theoretical separatrix of the resonance in the resonantintegral approximation (in reality, instead of the separatrix we have a thin separatrix-likechaotic layer). For the initial conditions taken along the line B ( u − u ), the orbit for whichΨ is maximum corresponds to a level curve tangent to the line, labeled Ψ ∗ . The initialcondition for u along this line, u , represents a good approximation to the exact resonantposition u res . The two values of u on the line B ( u − u ) satisfying Ψ = Ψ mnx correspond tothe intersection of the separatrix with the line B ( u − u ) (∆ u min and ∆ u min , in blue), andprovide an estimation of the width of the resonance. Right panel - Values of the resonantintegral Ψ along the line B ( u − u ). The position of the maximum of the function correspondsto u (green dot). The value of Ψ mnx (black line) defines the position of the two borders ofthe resonance ∆ u min and ∆ u max (blue dots). On the other hand, the points of intersection of the line Ψ = Ψ mnx with the curve of the resonantintegral mark the values u , u , and hence ∆ u min = u − u , ∆ u max = u − u , where the theoreticalseparatrix intersects the plane of the stability map. The corresponding values of e p can be foundthrough e p ,i = (cid:104) − Y (cid:16) u i , v i = B ( u i − u ) , φ f ; E (cid:17)(cid:105) / , with i = 1 , E allows to obtain the wholelocus of the theoretical center as well as the theoretical boundary of the secondary resonance onthe FLI stability map. Figure 3 shows an example of the location of the center and borders of asecondary resonance, with the method of the resonant integral, for the case of the 1:6 secondaryresonance ( µ = 0 . e (cid:48) = 0 .
02. The position of the center of the resonance is denoted by adashed line, and the inner and outer borders are denoted by thick solid lines. By comparison with theunderlaying FLI stability map, we can see that both the center of the resonance and the outer border∆ u max are understimated by this computation, proving that the overall estimation of the resonancewidth is not accurate. On the other hand, the key remark is that the method turns to be extremelyefficient in the location of the inner border. The approximate position of ∆ u min is well determinedin the whole range of proper eccentricity values considered 0 < e p, < . µ = 0 . µ = 0 . µ = 0 . µ = 0 . e (cid:48) = 0 .
02. In all the panels, the location of the inner border ∆ u min is µ = 0 . B = 0 and e (cid:48) = 0 .
02. The solid lines correspond to the inner an outer border ofthe resonance, the dashed line correspond to the estimation of the center of the resonance.The underlying image gives the numerical stability map, using the FLI value in grayscale. shown with a thick black line on top of the corresponding FLI stability map. We observe that thislimit divides the space of proper elements in two regions: the inner domain from ∆ u = 0 to ∆ u min is populated mainly by regular orbits, and exhibits also some isolated resonances of small width,in which the orbits can only be weakly chaotic and remain practically stable. On the contrary, thedomain external to ∆ u min is dominated by the presence of conspicuous resonances as well as regions ofstrong chaos. It is remarkable that the analytical determination of the inner border of the resonances,which is based on an integrable approximation to the Hamiltonian (i.e. the resonant normal form),can still provide an accurate limit even in domains of the phase space where the resonant orbits are,in reality, chaotic. It is this robustness of the inner border determination which renders the wholeapproach useful in practice.Figures 5 and 6 show, now, more examples of the applicability as well as the level of approximationof the method. Figure 5 shows the stability maps for µ = 0 . e (cid:48) = 0 .
02 (panel a), e (cid:48) = 0 .
06 (panel b), e (cid:48) = 0 . e (cid:48) = 0 (dotted thin line) and e (cid:48) = 0 .
02 (thick line) in (a), e (cid:48) = 0 (dotted thin line) and e (cid:48) = 0 .
06 (thick line) in (b) and e (cid:48) = 0 (dotted thin line) and e (cid:48) = 0 . e (cid:48) = 0 to only e (cid:48) = 0 . u (cid:39) .
4. In fact,we observe that, with increasing e (cid:48) , so called ‘transverse’ resonances, i.e. involving also the secularfrequency g , i.e. of the form m f ω f + m s ω s + m g g = 0 with m g (cid:54) = 0, appear near this border. Forexample, the 1:6:1 resonance at u ∼ .
25 in panel (a) of Fig. 5 moves towards the border at u ≈ . e (cid:48) these resonances have a small width and remain isolated within the inner stability domain, while, as e (cid:48) increases, all resonances (main or transverse) grow in size and move outwards, until they enter tothe region of strong chaos. As revealed in the panels of Fig. 5, these two effects (the moving of theresonances outwards and the refilling of the stable region with transverse resonances) counteract eachother in such a way that the border separating the inner domain of stability from the outer chaotic u min , thick black line), for thesecondary resonances 1:5 ( µ = 0 . B = 0 .
03, panel a), 1:6 ( µ = 0 . B = 0, panel b),1:7 ( µ = 0 . B = 0 . µ = 0 . B = 0, panel d), and e (cid:48) = 0 . domain remains practically in the same place. Due to this effect, we can see that even the estimationof the border via the resonant normal form corresponding to the circular case ( e (cid:48) = 0, dotted thinline) practically suffices to obtain a good approximation of the border of the effective stability domain.Also, regarding the Trojan’s body eccentricity, parameterized by e p, , one remarks that stable domainsof all the secondary resonances, beyond the main stability domain, survive only for small values of e p, .This is because the amplitude of the separatrix pulsation increases as the eccentricity of the Trojanbody increases. As a consequence, we find that the border of the main domain of stability is moresharp, and, thus, in general, better represented by the analytical resonance limit as e p, increases.Similar results are found in Fig. 6, showing the stability maps for µ = 0 . e (cid:48) = 0 .
02 (panel a), e (cid:48) = 0 .
08 (panel b). The estimated borders are found by the resonant normal form determination for µ = 0 . e (cid:48) = 0 (circular case, dotted thin line) and e (cid:48) = 0 .
02 (thick line)in (a), and e (cid:48) = 0 (dotted thin line) and e = 0 .
08 (thick line) in (b). The margin between the twotheoretical curves is again small (of about 0 .
02 rad in ∆ u ), while, again, the determination of theborder of the stability domain using the circular model suffices to practically obtain an accurate limitof the domain of stability. In fact, in both Figures 5 and 6 the extent occupied by the stable partsof the corresponding resonances is determined by the separatrix pulsation effect. The amplitude ofthe pulsation depends on terms absent from the ‘basic model’, thus this effect cannot be modelledusing only the resonant integrals of the basic model. However, as a rule of thumb we find that theborder of the domain of stability lies always between two theoretical border determinations by thebasic model, i.e., one using the circular model e (cid:48) = 0 and a second using a moderate value of the µ = 0 . B = 0) and threevalues of the eccentricity e (cid:48) = 0 .
02 (a), e (cid:48) = 0 .
06 (b), e (cid:48) = 0 . u min for the parameters µ = 0 . e (cid:48) = 0 (circular case) in all three panels, while the thick line corresponds to e (cid:48) = 0 .
02 in (a), e (cid:48) = 0 .
06 in (b) and e (cid:48) = 0 . primary’s eccentricity, e.g. e (cid:48) = 0 . The investigation in the previous subsection focused on particular values of µ selected with the criterionthat, for low eccentricities of either the primary perturber or the test body ( e (cid:48) , e p, < . n with n = 5 , , ... . Repeating a comparison between FLI maps and innermost separatrix bordersof secondary resonances, a behavior similar to Figures 6 (resonance 1:5, for µ = 0 . µ = 0 . µ = 0 . µ = 0 . µ = 0 . µ = 0 . µ = 0 . µ = 0 . µ are shifted positively with respect to the bifurcation values µ = µ n of each corresponding 1: n short period family in the basic model. The shift reflects the fact that,keeping e (cid:48) , e p, constant, and increasing µ as µ = µ n + ∆ µ , with ∆ µ >
0, the resonance 1: n movesoutwards from the libration center, i.e., towards higher libration amplitudes ∆ u , as ∆ µ increases.In the resonant integral approximation, the outward motion of each resonance is accompanied by anincrease of its separatrix width. However, the integrable aproximation fails due to resonance overlapwith nearby resonances as ∆ µ increases. This antagonism between outward expansion and resonanceoverlap determines the real limit of the domain of stability (see [40], [6] for a description of thisphenomenon in simple dynamical maps).The bifurcation value µ m : n for the m : n short-period family of the basic model can be estimatedby the root for µ of the equation: m (1 − µ/
8) = n (cid:115) µ (cid:18)
98 + 63 e (cid:48)
16 + 129 e p (cid:19) (36)Applying Eq. (36) to the 1:6 resonance, we find µ ≈ . e (cid:48) = e p = 0 .
02, while µ ≈ . e (cid:48) = e p = 0 .
1. As evident from Fig. 5, the resonance is clearly dominant at µ = 0 . µ ≥ . µ = 0 . µ ≈ . µ ≈ . e (cid:48) = e p, = 0 .
02, while µ ≈ . e (cid:48) = e p, = 0 .
1. Thus µ − µ ≈ . µ µ = 0 . B = 0 .
03) andtwo values of the eccentricity e (cid:48) = 0 .
02 (a), e (cid:48) = 0 .
08 (b). The dotted thin line correspondsto the analytical determination of ∆ u min for the parameters µ = 0 . e (cid:48) = 0 (circularcase) in the two panels, while the thick line corresponds to e (cid:48) = 0 .
02 in (a) and e (cid:48) = 0 .
08 in(b). separating the resonances 1:6 and 1:5 is about 3-4 times larger than the interval of values of ∆ µ forwhich the validity of the resonant integral computation using one particular resonance is satisfactory.In principle, in order to bridge the gap between the two resonances, one has to use higher orderresonances of the basic model, since the border of the domain of stability is always delimited by onesuch resonance. In practice, we find that it suffices to consider the basic resonances 1: n and their firstFarey tree combination, i.e., the resonances 2:(2 n −
1) which bifurcate at intermediate values of µ , i.e. µ n < µ n − < µ n − for fixed e (cid:48) , e p . Figure 7 exemplifies the transition from the dominance ofthe 1:6 to the 1:5 resonance via the 2:11 resonance of the basic model, for two values of the primary’seccentricity e (cid:48) = 0 .
02 (upper row) and e (cid:48) = 0 .
08 (lower row). The collapse of the inner bordercalculation for the 1:6 resonances starts near µ = 0 . µ = 0 . ≈ µ = 0 . µ the 1:5 secondary resonance of the basic model bifurcatesfor large enough values of the eccentricities, a fact which implies that the whole domain beyond theinnermost separatrix of the 1:5 resonance should now be considered as outside the main stabilitydomain. Indeed, although these secondary resonances are still very stable for very low eccentricities,we see that they essentially disappear for values of the eccentricities near ≈ . µ ( µ = 0 . transverse in [33], i.e. resonances involving all three short, synodic and secular frequencies. In particular, we seethat the border of stability, which for µ = 0 . µ increases, for higher values of the eccentricities. This effect leaves small windowsof values of µ for which, for low eccentricities, the border of stability may appear dominated bythe innermost separatrix of some transverse resonance (e.g. the resonance 1:5:1 in panel (c) for µ = 0 . e (cid:48) = 0 . µ , the innermost separatrix border of the1:5 resonance appears also in the upper part of the stability map for higher primary’s eccentricity,i.e., e (cid:48) = 0 .
08 (panel g). As a consequence, although a clear dominance of the 1:5 resonance occursfor all eccentricities beyond µ = 0 . µ = 0 . e (cid:48) = 0 .
02 (a), µ = 0 . e (cid:48) = 0 .
02 (b), µ = 0 . e (cid:48) = 0 .
02 (c), µ = 0 . e (cid:48) = 0 .
02 (d), µ = 0 . e (cid:48) = 0 .
02 (e), µ = 0 . e (cid:48) = 0 .
08 (f), µ = 0 . e (cid:48) = 0 .
08 (g), µ = 0 . e (cid:48) = 0 .
08 (h), µ = 0 . e (cid:48) = 0 .
08 (i), and µ = 0 . e (cid:48) = 0 .
08 (j). The thick lines yield theanalytical estimation of the border of stability for the corresponding value of µ and e (cid:48) = 0 . e (cid:48) = 0 .
08 for the lower row panels. The dotted thin line yieldsthe analytical estimation of the border for the corresponding value of µ and e (cid:48) = 0 (circularcase) in all the panels. B = 0 .
03 for all the panels. wide range of eccentricities already at µ = 0 . n or 2:(2 n − n integer. As a rule of thumb, for given parameters µ, e (cid:48) we choose the limiting resonance as the rational number closer to the frequency ratio: f = (cid:114) µ (cid:16) + e (cid:48) + e p (cid:17) (1 − µ/
8) (37)for values of e p within our domain of interest. As an additional test, we examine the robustness of the above results against changing the dynamicalmodel for Trojan orbits. Several formation scenaria discussed in literature ([1], [3], [16], [25], [32])have allowed relatively massive Trojan planets (of mass ∼ µ varies from µ = 0 . µ = 0 . µ = 0 . e (cid:48) = 0 .
02, (b). µ = 0 . e (cid:48) = 0 .
02, (c). µ = 0 . e (cid:48) = 0 .
02, (d). µ = 0 . e (cid:48) = 0 . µ = 0 . e (cid:48) = 0 .
08, (f). µ = 0 . e (cid:48) = 0 .
08, (g). µ = 0 . e (cid:48) = 0 .
08, H. µ = 0 . e (cid:48) = 0 .
08. The domains of various resonances (including transverse ones) aremarked in the same plots.
As an example, Figure 9 compares the border of stability in the planar ERTBP for µ = 0 . H = p m + p m + ( p + p ) m − Gm m r − Gm m r − Gm m ∆ (38)where m , m , m are the masses of the star, perturbing primary and Trojan planet respectively, r , r the heliocentric positions of the primary and Trojan planet, ∆ = r − r and p , b their correspondingbarycentric momenta. In order to use same units as in the ERTBP, one notes that the equations ofmotion in Cartesian heliocentric co-ordinates r ≡ ( r x , r y ), r ≡ ( r x , r y ), and barycentric velocities p /m ≡ ( v x , v y ), p /m ≡ ( v x , v y ) only depend on the variables r ix , r iy , v ix , v iy , i = 1 , Gm , Gm , Gm . Then, we solve Gm + Gm = 1, Gm µ and assign a value to Gm = Gm ( m /m ) according to the considered mass ratio µ = m /m .In order now to obtain comparable FLI maps in the two problems, we proceed as follows. Forevery point on the plane of the stability map of the ERTBP (such as in Fig. 9a), we compute thecorresponding heliocentric positions and velocities of both the primary and the Trojan, i.e. ( r ix ( t = 0), r iy ( t = 0), ˙ r ix ( t = 0) and ˙ r iy ( t = 0), with i = 1 ,
2. From Hamilton’s equations of (38) one readilysees that the barycentric velocities ( v ix , v iy ), i = 1 , r ix ( t = 0), r iy ( t = 0), v ix ( t = 0) and v iy ( t = 0) needed in orderto integrate the full Three Body problem. Via the same process we assign also corresponding initialconditions for the variational equations of motion in the two problems.Figure 9 shows the comparison of the FLI stability maps in the case of dominance of the 1:6resonance, at µ = 0 . µ evolves, i.e., µ = 3 × − (1 Earth mass, upper row), or µ = 3 × − (10 Earth masses, lower row). The left, middle and right panels correspond to initial eccentricitiesof the primary equal to e (cid:48) = 0 . e (cid:48) = 0 .
06 and e (cid:48) = 0 .
1. Thus, these maps are comparable withthe ones under the ERTBP (Fig. 5). The main observation is that switching on the mass m resultsin a considerable reduction of the area occupied by the stable domains of the secondary resonances. µ (cid:54) = 0assigned to the Trojan body. The mass of the primary is always µ = 0 . e (cid:48) of the primary and mass µ of the Trojan) are(a). e (cid:48) = 0 . µ = 3 × − , (b). e (cid:48) = 0 . µ = 3 × − , (c). e (cid:48) = 0 . µ = 3 × − , (d). e (cid:48) = 0 . µ = 3 × − , (e). e (cid:48) = 0 . µ = 3 × − , (f). e (cid:48) = 0 . µ = 3 × − . Theanalytical curves are those of Fig. 5. This is mostly caused by the secular variations induced on the orbit of the primary planet, whichincrease the amplitude of modulation of the separatrices of each secondary resonance. However, themain domain of stability remains nearly unaffected by these phenomena, and retains a quite similarwidth in all simulations with different masses µ . We only see some transverse resonances penetratingthe lowermost (with respect to the eccentricities) part of the stability map for µ as large as 10 Earthmasses. On the other hand, the analytical determination of the innermost separatrix via the resonantintegral of the ‘basic model’ yields an estimate of the border of the main stability domain whichremains robust against the increase of µ . In the present work, we discussed a new application for the ‘basic Hamiltonian model’ H b for Trojanmotions presented originally in [33]: this is the determination of the border of effective stability, usingthe theoretical separatrices of the most conspicuous secondary resonances of H b . In detail:1) We compute resonant normal forms for various secondary resonances of H b , using an ‘asym-metric expansion’ for the Hamiltonian (see Section 2), which allows to speed up the convergenceof both the original polynomial representation of the Hamiltonian as well as its normal form. Theimprovement obtained by the asymmetric expansion is demonstrated with numerical examples.2) Using the classical normal form construction with Lie series in order to compute a resonantnormal form for a specific secondary resonance, one ends with an expression for an invariant of thenormal form called the ‘resonant integral’ Ψ (see Section 3). The level curves of Ψ allow, in turn, toobtain a theoretical phase portrait on a surface of section, and in particular to compute theoreticalseparatrices as well as the center of the secondary resonance.3) The method typically yields underestimates of the position of the center and outer separatrix of the resonance, but a very efficient determination of the inner (closer to the libration center) separatrixof the resonance.4) We argued that the inner limit ∆ u min ( e p, ) found in this way represents a clear border whichexists in numerical stability maps between two well distinct domains in the space of the properelements (∆ u, e p, ) (see Section 4 for definitions). The inner domain is populated by regular orbitsand isolated resonances with regular or marginally chaotic orbits, while the outer domain hosts eitherclosely packed secondary resonances or a strongly chaotic domain. In fact, with increasing value of theprimary’s eccentricity e (cid:48) , a modulation mechanism essentially wipes out all the resonances, creating alarge outer domain of strong chaos. As a consequence, we argued that the inner domain, delimited bythe innermost theoretical separatrix of the most conspicuous secondary resonance of H b practicallycoincides with the limit of the effective stability domain for Trojan motions.5) We demonstrated that the role of the secondary resonances of the basic model H b , as delimitersof the domain of effective stability, covers most of the values of the parameters entering the problem(primary’s mass and eccentricity, Trojan body’s eccentricity), while it remains robust even in the fullThree Body problem, for Trojan bodies of mass ∼ Acknowledgements:
Useful discussions with Prof. U. Locatelli are gratefully acknowledged. R.I.P.was supported by the Research Comittee of the Academy of Athens, under the grant 200/854.
Appendix A
The variables corresponding to the three degrees of freedom appearing in the expression of the basicHamiltonian H b in Eq.(5), ( u, v ), ( Y f , φ f ) and ( Y p , φ p ) are given in terms of the orbital elements asfollows: u = λ − λ (cid:48) − π , (39) v = √ a − , (40) β = ω − φ (cid:48) ,y = √ a (cid:16)(cid:112) − e − (cid:17) ,V = (cid:112) − y sin β − (cid:112) − y sin β ,W = (cid:112) − y cos β − (cid:112) − y cos β ,Y = − (cid:18) W + V (cid:19) φ = arctan (cid:18) VW (cid:19) (41) φ f = λ (cid:48) − φ , (42) Y f = (cid:90) ∂E∂λ (cid:48) d t + v , (43) Y p = Y − Y f , (44)where λ , ω , a and e are the mean longitude, the longitude of the perihelion, the major semiaxis andeccentricity of the Trojan body, λ (cid:48) and φ (cid:48) = ω (cid:48) are the mean longitude and longitude of the perihelionof the perturber, β = π/ y = √ − e (cid:48) −
1, and E represents the total energy of the Trojan ascomputed from Eq. (1) (see [33] for further details in the construction). Appendix B
The asymmetric expansion in terms of u = τ − π/
3, up to a generic order K for the functions cos τ (2 − τ ) N/ , sin τ (2 − τ ) N/ , cos M τ and sin M τ , with N, M ∈ N fixed is given bycos τ (2 − τ ) N/ = 12 N/ K (cid:88) k =0 M ( k ) u k + O ( u K ) , where M ( k ) = K (cid:88) i = k i ! F ( i ) ( π/ (cid:18) ik (cid:19) (cid:16) − π (cid:17) i − k , sin τ (2 − τ ) N/ = 12 N/ K (cid:88) k =0 M ( k ) u k + O ( u K ) , where M ( k ) = K (cid:88) i = k i ! G ( i ) ( π/ (cid:18) ik (cid:19) (cid:16) − π (cid:17) i − k , cos M τ = K (cid:88) k =0 M ( k ) u k + O ( u K ) , where M ( k ) = K (cid:88) i = k i ! B ( i ) M,M (cid:18) ik (cid:19) (cid:16) − π (cid:17) i − k , sin M τ = K (cid:88) k =0 M ( k ) u k + O ( u K ) , where M ( k ) = K (cid:88) i = k i ! C ( i ) M,M (cid:18) ik (cid:19) (cid:16) − π (cid:17) i − k , and F ( n ) ( π/
2) = [ n − ] (cid:88) i =1 ( n, i −
1) ( − i f ( n − (2 i − ( π/ ,G ( n ) ( π/
2) = [ n ] (cid:88) i =0 ( n, i ) ( − i f ( n − i ) ( π/ , with [ n − ] the integer part of n − , and [ n ] the integer part of n ; the derivatives f ( n ) are given by f ( n ) ( π/
2) = n (cid:88) m =1 A ( n ) m,m ;the coefficients A ( n ) m,m , B ( n ) M,M and C ( n ) M,M are given by A ( n ) m,m = − A ( n − m,m − − (cid:18) m −
1) + N (cid:19) A ( n − m − ,m − , A (1)1 , = − N ,B ( n ) M,M = − B ( n − M,M − + ( M + 1) B ( n − M,M +1 , B (1)1 , = − M ,C ( n ) M,M = C ( n − M,M − − ( M + 1) C ( n − M,M +1 , C (1)1 , = M .
For a proof of these formulæ, we refer the reader to [36].
References [1] Beaug´e, C., S´andor, Z., ´Erdi, B., S¨uli, A., (2007) Co-orbital terrestrial planets in exoplanetarysystems: a formation scenario,
Astron. Astrophys. , p 359.[2] Chirikov, B.V., Lieberman, M.A., Shepelyansky, D.L., Vivaldi, F.M., (1985) A theory of modula-tional diffusion,
Physica D , p 289.[3] Cresswell, P., Nelson, R.P., (2009) On the growth and stability of Trojan planets, Astron. Astro-phys. , p 1141.[4] Dobrovolskis, A., (2013) Effects of Trojan exoplanets on the reflex motions of their parent stars,
Icarus , p 1635. [5] Dvorak, R., Bazs´o, A., Zhou, L.Y., (2010) Where are the Uranus Trojans?,
Celest. Mech. Dyn.Astron. , p 51.[6] Efthymiopoulos, C., Contopoulos, G., Voglis, N., (1999) Cantori, islands and asymptotic curvesin the stickiness region
Celest. Mech. Dyn. Astron. , p 221.[7] Efthymiopoulos, C., Giorgilli, A., Contopoulos, G., (2004) Nonconvergence on formal integrals:II. Improved estimates for the optimal order of truncations J. Phys. A , p 10831.[8] Efthymiopoulos, C. (2012) Canonical perturbation theory, stability and diffusion in Hamiltoniansystems: applications in dynamical astronomy, in Third La Plata International School on Astro-nomy and Geophysics: Chaos, diffusion and non-integrability in Hamiltonian Systems - Applica-tions to Astronomy , Cincotta, P.M., Giordano, C.M., Efthymiopoulos, C., eds, p 3.[9] Efthymiopoulos, C., (2013) High order normal form stability estimates for co-orbital motion,
Celest. Mech. Dyn. Astron. , p 101.[10] ´Erdi, B., (1988) Long periodic perturbations of Trojan asteroids,
Celest. Mech. Dyn. Astron. ,p 303.[11] ´Erdi, B., (1997) The Trojan Problem, Celest. Mech. Dyn. Astron. , p 149.[12] ´Erdi, B., Sandor, Z., (2005) Stability of co-orbital motion in exoplanetary systems, Celest. Mech.Dyn. Astron. , p 113.[13] ´Erdi, B., Nagy, I., S´andor, Z., S¨uli, A., Fr¨ohlich, G., (2007) Secondary resonances of co-orbitalmotions, MNRAS , p 33.[14] Froeschl´e, C., Guzzo, M., Lega, E., (2000) Graphical evolution of the Arnold web: from order tochaos,
Science , p 2108.[15] Giorgilli, A., Skokos, C., (1997) On the stability of the Trojan asteroids,
Astron. Astrophys. ,p 254.[16] Giuppone, C., Ben´ıtez-Llambay, P, Beaug´e, C., (2012) Origin and detectability of co-orbitalplanets from radial velocity data,
MNRAS , p 356.[17] Haghighipour, N., Capen, S., Hinse, T., (2013) Detection of Earth-mass and super-Earth Trojanplanets using transit timing variation method,
Celest. Mech. Dyn. Astron. , p 75.[18] Laughlin, G., Chambers, J.E., (2002) Extrasolar Trojans: the viability and detectability ofplanets in the 1:1 Resonance,
Astron.J. , p 592.[19] Leleu, A., Robutel, P., Correia, A.C.M. (2015) Detectability of quasi-circular co-orbital planets.Application to the radial velocity technique,
Astron.Astrophys. , p A128[20] Leleu, A., (2016) Dynamics of co-orbital exoplanets - Ph.D. Thesis, arXiv:1701.05585 [21] Leleu, A., Robutel, P., Correia, A.C.M., Lillo-Box, J. (2017) Detection of co-orbital planets bycombining transit and radial-velocity measurements,
Astron.Astrophys. , p L7[22] Levison, H., Shoemaker, E., Shoemaker, C., (1997) Dynamical evolution of Jupiter’s Trojanasteroids,
Nature , p 42.[23] Lhotka, C., Efthymiopoulos, C., Dvorak, R., (2008) Nekhoroshev stability at L L MNRAS , p 1165.[24] Lykawka, P.S., Horner, J., Jones, B.W., Mukai, T., (2011) Origin and dynamical evolution ofNeptune Trojans - II. Long term evolution,
MNRAS (1), p 537.[25] Lyra, W., Johansen, A., Klahr, H., Piskunov, N., (2009) Standing on the shoulders of giants:Trojan Earths and vortex trapping in low mass self-gravitating protoplanetary disks of gas andsolids,
MNRAS , p 1125.[26] Marzari, F., Scholl, H., (2007) Dynamics of Jupiter Trojans during the 2:1 mean motion resonancecrossing of Jupiter and Saturn,
MNRAS , p 479. [27] Milani, A., (1993) The Trojan asteroid belt: proper elements, stability, chaos and families,
Celest.Mech. Dyn. Astron. , p 59.[28] Morais, M.H.M, (1999) A secular theory for Trojan-type motion, Astron. Astrophys. , p 318.[29] Morais, M.H.M, (2001) Hamiltonian formulation on the secular theory for a Trojan-type motion,
Astron. Astrophys. , p 677.[30] Nauenberg, M., (2002) Stability and eccentricity for two planets in a 1:1 resonances, and theirpossible occurrence in extrasolar planetary systems,
Astron. J. , p 2332.[31] Neishtadt, A.I., (1987) On the change in the adiabatic invariant on crossing a separatrix insystems with two degrees of freedom,
Prikl. Matem. Mekhan. (5), p. 750; PMM USSR (5), p586.[32] Pierens, A., Raymond, S.N., (2014) Disruption of co-orbital (1:1) planetary resonances duringgas-driven orbital migration, MNRAS , (2), p 2296.[33] P´aez, R.I., Efthymiopoulos, C., (2015) Trojan resonant dynamics, stability, and chaotic diffusion,for parameters relevant to exoplanetary systems, Celest. Mech. Dyn. Astron. , (2), p 139.[34] P´aez, R.I., Locatelli, U., (2015) Trojan dynamics well approximated by a new Hamiltoniannormal form, MNRAS (2), p 2177.[35] P´aez, R.I., Locatelli, U., Efthymiopoulos, C. (2016)
Celest.Mech.Dyn.Astron , p 519[36] P´aez, R.I. (2016) New normal form approaches adapted to the Trojan problem - Ph.D. Thesis, arXiv:1703.08819 [37] Robutel, P., Gabern, F., (2006) The resonant structure of Jupiter’s Trojan asteroids - I. Longterm stability and diffusion,
MNRAS , p 1463.[38] Schwarz, R., S¨uli, ´A., Dvorak, R., Pilat-Lohinger, E., (2009) Stability of Trojan planets inmultiplanetary systems,
Celest. Mech. Dyn. Astron. , p 69.[39] Tsiganis, K., Varvoglis, H., Dvorak, R., (2005) Chaotic diffusion and effective stability of JupiterTrojans,
Celest. Mech. Dyn. Astron. , p 71.[40] Voglis, N., Efthymiopoulos, C., (1998) Angular dynamical spectra. A new method for determiningfrequencies, weak chaos and cantori J. Phys. A , p. 2913., p. 2913.