On the extension of the Laplace-Lagrange secular theory to order two in the masses for extrasolar systems
aa r X i v : . [ m a t h - ph ] J un Celestial Mechanics and Dynamical Astronomy
DOI 10.1007/s10569-013-9501-z
On the extension of the Laplace-Lagrange secular theoryto order two in the masses for extrasolar systems
Anne-Sophie Libert · Marco Sansottera
Received: 9 November 2012 / Accepted: 12 June 2013
Abstract
We study the secular evolution of several exoplanetary systems by extend-ing the Laplace-Lagrange theory to order two in the masses. Using an expansion ofthe Hamiltonian in the Poincar´e canonical variables, we determine the fundamentalfrequencies of the motion and compute analytically the long-term evolution of the ke-plerian elements. Our study clearly shows that, for systems close to a mean-motionresonance, the second order approximation describes their secular evolution more ac-curately than the usually adopted first order one. Moreover, this approach takes intoaccount the influence of the mean anomalies on the secular dynamics. Finally, we setup a simple criterion that is useful to discriminate between three different categories ofplanetary systems: (i) secular systems (HD 11964, HD 74156, HD 134987, HD 163607,HD 12661 and HD 147018); (ii) systems near a mean-motion resonance (HD 11506,HD 177830, HD 9446, HD 169830 and υ Andromedae); (iii) systems really close to or in a mean-motion resonance (HD 108874, HD 128311 and HD 183263). Keywords extrasolar systems · n-body problem · secular dynamics · normal formsmethod · proximity to mean-motion resonances The study of the dynamics of planetary systems is a long standing and challengingproblem. The classical perturbation theory, mainly developed by Lagrange and Laplace,uses the circular approximation as a reference for the orbits. The discovery of extrasolar
A.-S. LibertnaXys, Department of Mathematics, University of Namur, 8 Rempart de la Vierge,5000 Namur, BelgiumObservatoire de Lille (LAL-IMCCE), CNRS-UMR8028, 1 Impasse de l’Observatoire,59000 Lille, FranceE-mail: [email protected]. SansotteranaXys, Department of Mathematics, University of Namur, 8 Rempart de la Vierge,5000 Namur, BelgiumE-mail: [email protected] planetary systems has opened a new field in Celestial Mechanics and nowadays morethan 100 multi-planetary systems have been discovered. In contrast to the Solar system,where the orbits of the planets are almost circular, the exoplanets usually describetrue ellipses with high eccentricities. Thus, the applicability of the classical approach,using the circular approximation as a reference, can be doubtful for these systems. Inthis work we revisit the classical Laplace-Lagrange theory for the secular motions ofthe pericenters of the planetary orbits, based only on a linear approximation of thedynamical equations, by considering also higher order terms.Previous works of Libert & Henrard (2005, 2006) for coplanar systems have gen-eralized the classical expansion of the perturbative potential to a higher order in theeccentricities, showing that this analytical model gives an accurate description of thebehavior of planetary systems which are not close to a mean-motion resonance, up tosurprisingly high eccentricities. Moreover, they have shown that an expansion up toorder 12 in the eccentricities is usually required for reproducing the secular behaviorof extrasolar planetary systems. This expansion has also been used by Beaug´e et al.(2006) to successfully reproduce the motions of irregular satellites with eccentricitiesup to 0 .
7. Veras & Armitage (2007) have highlighted the limitations of lower orderexpansions; using only a fourth-order expansion in the eccentricities, they did not re-produce, even qualitatively, the secular dynamics of extrasolar planetary systems. Allthe previous results have been obtained considering a secular Hamiltonian at orderone in the masses. Let us also remark that an alternative octupole-level secular theoryhas been developed for systems that exhibit hierarchical behavior (see, e.g., Ford et al.(2000), Lee & Peale (2003), Naoz et al. (2011), Katz et al. (2011) and Libert & Delsate(2012)). However, this approach is not suitable for our study, as we will also considersystems with large semi-major axes ratio.Considering the secular dynamics of the Solar system, Lagrange and Laplace showedthat the proximity to the 5:2 mean-motion resonance between Jupiter and Saturn leadsto large perturbations in the secular motion, explaining the so-called “great inequal-ity”. Let us stress that when referring to secular evolution we mean long-term evo-lution, possibly including the long-term effects of near resonances. Indeed the termsof the perturbation associated to mean-motion resonances have small frequencies andthus influence the secular behavior of the system. A good description of the seculardynamics of an exoplanetary system should include the effects of the mean-motionresonances on their secular long-term evolution. Therefore, we replace the classical cir-cular approximation with a torus which is invariant up to order two in the masses,this is the so-called Hamiltonian at order two in the masses. The benefit of a secondorder approach has been clearly highlighted in Laskar (1988), see Table 2 in that paper,where a comparison between the fundamental frequencies of the planetary motion ofthe Solar system, is given for different approximations. Concerning the problem of thestability of the Solar system, the celebrated theorems of Kolmogorov and Nekhoroshevallowed to make substantial progress. Nevertheless, in order to apply these theorems,a crucial point is to consider the secular Hamiltonian at order two in the massesin order to have a good approximation of the secular dynamics (e.g., this allows, inLocatelli & Giorgilli (2007), to deal with the true values of the planetary masses, whilethe first order approximation used in Locatelli & Giorgilli (2005) forces the authorsto reduce the masses of the planets by a factor 10). In recent years, the estimates forthe applicability of both Kolmogorov and Nekhoroshev theorems to realistic modelsof some part of the Solar system have been improved by some authors (e.g., Robutel (1995), Fejoz (2005), Celletti & Chierchia (2005), Gabern (2005), Locatelli & Giorgilli(2007), Giorgilli et al. (2009) and Sansottera et al. (2013, 2011)).In the present contribution, we study the secular dynamics of extrasolar planetarysystems consisting of two coplanar planets. The aim is to reconstruct the evolutionof the eccentricities and pericenters of the planets by using analytical techniques, ex-tending the Laplace-Lagrange theory to order two in the masses. To do so, we extendthe results in Libert & Henrard (2005, 2006), replacing the first order averaged Hamil-tonian, with the one at order two in the masses, and show the improvements of thisapproximation on the study of the secular evolution of extrasolar systems. In particular,we determine the fundamental frequencies of the motion and compute analytically thelong-term evolution of the keplerian elements. Furthermore, we show that the Hamil-tonian at order two in the masses describes accurately the secular dynamics of systemsclose to a mean-motion resonance and, as a byproduct, we also give an estimate of theproximity to a mean-motion resonance of the two-planet extrasolar systems discoveredso far.The paper is organized as follows. In Section 2, we describe the expansion of theHamiltonian of the planar three-body problem in Poincar´e variables. Following theLagrange approach, we focus, in Section 3, on the secular part of the Hamiltonian andderive the secular Hamiltonian at order two in the masses. In Section 4, we constructa high-order Birkhoff normal form, using the Lie series method, that leads to a verysimple form of the equations of motion, being function of the actions only. We also showhow to compute the secular frequencies and perform a long-term analytical integrationof the motion of the planets. In Section 5, we apply our method to the υ Andromedaesystem and show that the second order approximation is well suited for systems closeto a mean-motion resonance. Furthermore, the influence of the mean anomaly on thelong-term evolution is pointed out in Section 6. In Section 7, we set up a criterion toevaluate the proximity of planetary systems to mean-motion resonances, and apply it tothe two-planet extrasolar systems discovered so far. Finally, our results are summarizedin Section 8. An appendix containing the expansion of the secular Hamiltonian of the υ Andromedae extrasolar system, up to order 6, follows.
We consider a system of three coplanar point bodies, mutually interacting according toNewton’s gravitational law, consisting of a central star P of mass m and two planets P and P of mass m and m , respectively. The indices 1 and 2 refer to the inner andouter planet, respectively.Let us now recall how the classical Poincar´e variables can be introduced to performthe expansion of the Hamiltonian around circular orbits. We follow the formalism intro-duced by Poincar´e (see Poincar´e (1892, 1905); for a modern exposition, see, e.g., Laskar(1989) and Laskar & Robutel (1995)). To remove the motion of the center of mass, weadopt the heliocentric coordinates † , r j = −−−→ P P j , with j = 1 ,
2. Denoting by ˜ r j themomenta conjugated to r j , the Hamiltonian of the system has four degrees of freedom,and reads F ( r , ˜ r ) = T (0) (˜ r ) + U (0) ( r ) + T (1) (˜ r ) + U (1) ( r ) , (1) † Let us note that the Jacobi variables are less suitable for our purpose, as they require aTaylor expansion in the planetary masses. where T (0) (˜ r ) = 12 X j =1 k ˜ r j k (cid:18) m + 1 m j (cid:19) , T (1) (˜ r ) = ˜ r · ˜ r m ,U (0) ( r ) = −G X j =1 m m j k r j k , U (1) ( r ) = −G m m k r − r k . The plane set of Poincar´e canonical variables is introduced as Λ j = m m j m + m j q G ( m + m j ) a j , λ j = M j + ω j ,ξ j = q Λ j r − q − e j cos ω j , η j = − q Λ j r − q − e j sin ω j , (2)for j = 1 , a j , e j , M j and ω j are the semi-major axis, the eccentricity, themean anomaly and the longitude of the pericenter of the j -th planet, respectively. Oneimmediately sees that both ξ j and η j are of the same order as the eccentricity e j . Usingthe Poincar´e canonical variables, the Hamiltonian becomes F ( Λ , λ , ξ , η ) = F (0) ( Λ ) + F (1) ( Λ , λ , ξ , η ) , (3)where F (0) = T (0) + U (0) is the keplerian part and F (1) = T (1) + U (1) the perturbation.Let us emphasize that the ratio F (1) /F (0) = O ( µ ) with µ = max { m /m , m /m } .Therefore, the time derivative of each variable is of order µ , except for λ . For thisreason we will refer to ( Λ , λ ) as the fast variables and to ( ξ , η ) as the secular variables .We proceed now by expanding the Hamiltonian (3) in Taylor-Fourier series. Wepick a fixed value Λ ∗ of the fast actions † and perform a translation, T F , defined as L = Λ − Λ ∗ . This canonical transformation leaves the coordinates λ , ξ and η unchanged. The trans-formed Hamiltonian H ( T ) = T F ( F ) can be expanded in power series of L , ξ and η around the origin. Forgetting an unessential constant, we rearrange the Hamiltonianof the system as H ( T ) = n ∗ · L + ∞ X j =2 h (Kep) j , ( L ) + µ ∞ X j =0 ∞ X j =0 h ( T ) j ,j ( L , λ , ξ , η ) , (4)where the functions h ( T ) j ,j are homogeneous polynomials of degree j in the fast actions L , degree j in the secular variables ( ξ , η ), and depend analytically and periodically onthe angles λ . The terms h (Kep) j , of the keplerian part are homogeneous polynomials ofdegree j in the fast actions L . We also expand h ( T ) j ,j in Fourier series of the angles λ .In the latter equation, the term which is both linear in the actions and independent ofall the other canonical variables (i.e., n ∗ · L ) has been isolated in view of its relevancein perturbation theory, as it will be discussed in the next section. † We recall that, as shown by Poisson, the semi-major axes are constant up to the secondorder in the masses. Here we expand around their initial values, but we could also have takentheir average values over a long-term numerical integration (see, e.g., Sansottera et al. (2013)).
All the expansions were carried out using a specially devised algebraic manipulator(see Giorgilli & Sansottera (2011)). In our computations we truncate the expansion asfollows. The keplerian part is expanded up to the quadratic terms. The terms h ( T ) j ,j include the linear terms in the fast actions L , all terms up to degree 12 in the secularvariables ( ξ , η ) and all terms up to the trigonometric degree 12 with respect to theangles λ . The choice of the limits in the expansion is uniform for all the systemsthat will be considered. However, as we will explain in the next section, the actuallimits for the computation of the secular approximation will be chosen as the lowestpossible orders in λ and ( ξ , η ), so as to include the main effects of the proximity to amean-motion resonance. In this section we discuss the procedure for computing the secular Hamiltonian viaelimination of the fast angles. The classical approach, usually found in the literature,consists in replacing the Hamiltonian H ( T ) , defined in (4), by H = 14 π Z π Z π H ( T ) d λ d λ . (5)The idea is that the effects due to the fast angles are negligible on the long-termevolution and this averaged Hamiltonian represents a “good approximation” of thesecular dynamics. This approach has been critically considered by Arnold, quotinghis book (i.e., Arnold (1989), Chapter 10) “this principle is neither a theorem, anaxiom, nor a definition, but rather a physical proposition, i.e., a vaguely formulatedand, strictly speaking, untrue assertion. Such assertion are often fruitful sources ofmathematical theorems.” .The secular Hamiltonian obtained in this way is the so-called approximation at order one in the masses (or “averaging by scissors”) and is the basis of the Laplace-Lagrange theory for the secular motion of perihelia and nodes of the planetary orbits.This approximation corresponds to fixing the value of Λ , that remains constant underthe flow, and thereby the semi-major axes. The averaged Hamiltonian, depending onlyon the secular variables, reduces the problem to a system with two degrees of freedom.Let us remark that the Laplace-Lagrange secular theory was developed just con-sidering the linear approximation of the dynamical equations. An extension of theLaplace-Lagrange theory for extrasolar systems, including also terms of higher orderin the eccentricities, can be found in Libert & Henrard (2005, 2006), where the au-thors show that a secular Hamiltonian at order one in the masses gives an accuratedescription of the long-term behavior for systems which are not close to a mean-motionresonance.One of the main achievements of the Laplace-Lagrange secular theory is the expla-nation of the “great inequality” between Jupiter and Saturn. Indeed, they have shownthat the near commensurability of the two mean-motions (the 5:2 near resonance) hasa great impact on the long-term behavior of the Solar system. For that reason, a gooddescription of the secular dynamics of an exoplanetary system should include a carefultreatment of the influence of mean-motion resonances on the long-term evolution.Our purpose is to consider a secular Hamiltonian at order two in the masses . Theidea is to remove perturbatively the dependency on the fast angles from the Hamilto-nian (4), considering terms up to order two in the masses. This can be done using the classical generating functions of the Hamilton-Jacobi formalism. Here instead we usethe Lie series formalism and implement the procedure in a way that takes into accountthe Kolmogorov algorithm for the construction of an invariant torus (see, Kolmogorov(1954)). This is only a technical point and does not affect the results, our choice is aquestion of convenience since the Lie series approach is a direct method and is muchmore effective from the computational point of view (see, e.g., Giorgilli & Locatelli(2003)).3.1 Approximation at order two in the massesLet us recall that in the expansion of the Hamiltonian H ( T ) , see (4), the perturbation isof order one in the masses and it is a polynomial in L , ξ and η , and a trigonometric poly-nomial in the fast angles λ . We remove part of the dependence on the fast angles per-forming a “Kolmogorov-like” normalization step, in the following sense. The suggestionof Kolmogorov is to give the Hamiltonian the normal form H ( L , λ ) = n ∗ · L + O ( L ),for which the existence of an invariant torus L = 0 is evident (where we consider thesecular variables just as parameters). We give the Hamiltonian the latter form up toterms of order O ( µ ). More precisely, we perform a canonical transformation whichremoves the dependence on the fast angles from terms which are independent of andlinear in the fast actions L (i.e., equations (6) and (7), respectively). Therefore wereplace the circular orbits of the Laplace-Lagrange theory with an approximate in-variant torus, thus establishing a better approximation as the starting point of theclassical theory. We also take into account the effects of near mean-motion resonancesby including in the averaging process the corresponding resonant harmonics, as willbe detailed hereafter. The procedure is a little cumbersome, and here we give only asketch of the main path. For a detailed exposition one can refer to Locatelli & Giorgilli(2007) and Sansottera et al. (2013).The expansion of the Hamiltonian H ( T ) , see (4), in view of the d’Alembert rules(see, e.g., Poincar´e (1905); see also Kholshevnikov (1997, 2001) for a modern approach),contains only specific combinations of terms. Let us consider the harmonic k · λ , where k = ( k , k ), and introduce the so-called “characteristic of the inequality” C I ( k ) = k + k . The degree in the secular variables of the non-zero terms appearing in the expansionmust have the same parity of C I ( k ) and is at least equal to |C I ( k ) | .It is well known that the terms of the expansion that have the main influence on thesecular evolution are the ones related to low order mean-motion resonances. Therefore,if the ratio n ∗ /n ∗ is close to the rational approximation k ∗ /k ∗ , then the effects due tothe harmonics ( k ∗ λ − k ∗ λ ) should be taken into account in the secular Hamiltonian.Let us also recall that the coefficients of the Fourier expansion decay exponentiallywith | k | = | k | + | k | , so we just need to consider low order resonances.Let us go into details. Consider a system close to the k ∗ : k ∗ mean-motion resonanceand define the vector k ∗ = ( k ∗ , − k ∗ ) and two integer parameters K F = | k ∗ | and K S = |C I ( k ∗ ) | . We denote by ⌈ f ⌉ λ ; K F the Fourier expansion of a function f truncatedin such a way that we keep only the harmonics satisfying the restriction 0 < | k | ≤ K F .The effect of the near mean-motion resonances is taken into account by choosing theparameters K S and K F as the lowest limits that include the corresponding resonantharmonics. Using the Lie series algorithm to calculate the canonical transformations (see, e.g., Henrard (1973) and Giorgilli (1995)), we transform the Hamiltonian (4) as b H ( O = exp L µ χ ( O H ( T ) , with the generating function µ χ ( O ( λ , ξ , η ) determinedby solving the equation X j =1 n ∗ j ∂ χ ( O ∂λ j + K S X j =0 l h ( T )0 ,j m λ ; K F ( λ , ξ , η ) = 0 . (6)Notice that, by definition, the average over the fast angles of l h ( T )0 ,j m λ ; K F is zero,which assures that (6) can be solved provided that the frequencies are non resonantup to order K F . The Hamiltonian b H ( O has the same form as H ( T ) in (4), with thefunctions h ( T ) j ,j replaced by new ones, ˆ h ( O j ,j , generated by expanding the Lie seriesexp L µ χ ( O H ( T ) and gathering all the terms having the same degree both in the fastactions and in the secular variables.We now perform a second canonical transformation H ( O = exp L µ χ ( O b H ( O ,where the generating function µ χ ( O ( L , λ , ξ , η ), which is linear with respect to L , isdetermined by solving the equation X j =1 n ∗ j ∂ χ ( O ∂λ j + K S X j =0 l ˆ h ( O ,j m λ ; K F ( L , λ , ξ , η ) = 0 . (7)Again, (7) can be solved provided the frequencies are non resonant up to order K F and the Hamiltonian H ( O can be written in a form similar to (4), namely H ( O ( L , λ , ξ , η ) = n ∗ · L + ∞ X j =2 h (Kep) j , ( L ) + µ ∞ X j =0 ∞ X j =0 h ( O j ,j ( L , λ , ξ , η ; µ ) + O ( µ ) , (8)where the new functions h ( O j ,j are computed as previously explained for ˆ h ( O j ,j and stillhave the same dependence on their arguments as h ( T ) j ,j in (4). As we are interested in asecond order approximation, we have neglected the contribution of the order O ( µ ) inthe canonical transformations associated to (6) and (7). Following a common practicein perturbation theory, we denote again by ( L , λ , ξ , η ) the transformed coordinates.In the following, we will denote by T O the canonical transformation induced bythe generating functions µ χ ( O and µ χ ( O , namely T O ( L , λ , ξ , η ) = exp L µ χ ( O ◦ exp L µ χ ( O ( L , λ , ξ , η ) . (9)Let us remark that the non resonant condition k · n ∗ = 0 , for 0 < | k | ≤ K F , does not imply that the canonical change of coordinates is convergent. Indeed, the terms k · n ∗ that appear as the denominators of the generating functions, even if they do notvanish, can produce the so-called small divisors . It is well known that the presence ofsmall divisors is a major problem in perturbation theory. Therefore, for each systemconsidered in this work, we check that the canonical transformation T O is near to theidentity and only in that case we proceed computing the approximation at order twoin the masses. H ( O in (8), we just need to perform an average overthe fast angles λ . More precisely, we consider the averaged Hamiltonian H (sec) ( ξ , η ) = D H ( O (cid:12)(cid:12) L =0 E λ , (10)namely we set L = 0 and average H ( O by removing all the Fourier harmonics de-pending on the angles. This results in replacing the orbit having zero eccentricity withan invariant torus of the unperturbed Hamiltonian. The Hamiltonian so constructed isthe secular one, describing the slow motion of the eccentricities and pericenters. Con-cerning the approximation at order one in the masses, let us recall that we directlyaverage the Hamiltonian H ( T ) , see equation (5).After the averaging over the fast angles, the secular Hamiltonian has two degreesof freedom and, in view of the d’Alembert rules, contains only terms of even degree in( ξ , η ). Therefore, the lowest order approximation of the secular Hamiltonian, namelyits quadratic part, is essentially the one considered in the Laplace-Lagrange theory.The origin ( ξ , η ) = (0 ,
0) is an elliptic equilibrium point, and it is well known thatone can find a linear canonical transformation ( ξ , η ) = D ( x , y ) which diagonalizes thequadratic part of the Hamiltonian, so that we may write H (sec) in the new coordinatesas H (sec) ( x , y ) = X j =1 ν j x j + y j H (0)2 ( x , y ) + H (0)4 ( x , y ) + . . . , (11)where ν j are the secular frequencies in the small oscillations limit and H (0)2 s is a homo-geneous polynomial of degree 2 s + 2 in ( x , y ) .To illustrate the transformations described hereabove, the secular Hamiltonian, H (sec) , of the υ Andromedae extrasolar system (see Section 5 for a detailed descriptionof the system and a discussion on its proximity to the 5:1 mean-motion resonance) isreported in appendix A. First and second order approximations in the masses, includingterms up to order 6 in ( ξ , η ), are given. Following Libert & Henrard (2006), we now aim to introduce an action-angle formu-lation, since its associated equations of motion are extremely simple and can be inte-grated analytically. The secular Hamiltonian (11) has the form of a perturbed systemof harmonic oscillators, and thus we can construct a Birkhoff normal form (see Birkhoff(1927)) introducing action-angle coordinates for the secular variables, by means of Lieseries (see, e.g., Hori (1966); Deprit (1969); Giorgilli (1995)). Finally, an analyticalintegration of the action-angle equations will allow us to check the accuracy of our sec-ular approximation, by comparing it to a direct numerical integration of the Newtonequations.4.1 Birkhoff normal form via Lie seriesAs the construction of the Birkhoff normal form via Lie series is explained in detail inmany previous studies, here we just briefly recall it, adapted to the present context.
First, we define a canonical transformation ( x , y ) = A ( I , ϕ ) introducing the usualaction-angle variables x j = p I j cos ϕ j , y j = p I j sin ϕ j , j = 1 , . (12)The secular Hamiltonian in these variables reads H (sec) ( I , ϕ ) = ν · I + H (0)2 ( I , ϕ ) + H (0)4 ( I , ϕ ) + . . . . (13)In order to remove the dependency of the secular angles ϕ in this expression, wecompute the Birkhoff normal form up to order r , H ( r ) = Z ( I ) + . . . + Z r ( I ) + R ( r ) ( I , ϕ ) , (14)where Z s , for s = 0 , . . . , r , is a homogeneous polynomial of degree s/ I and iszero for odd s . Only the remainder, R ( r ) ( I , ϕ ), depends also on the angles ϕ . Again,with a little abuse of notation, we denote by ( I , ϕ ) the new coordinates. At each order s >
0, we determine the generating function X ( s ) , by solving the equation n X ( s ) , ν · I o + H s ( I , ϕ ) = Z s ( I ) . (15)Using the Lie series, we calculate the new Hamiltonian as H ( s +1) = exp L X ( s +1) H ( s ) ,provided that the non-resonance condition k · ν = 0 for k ∈ Z such that 0 < | k | ≤ s + 2 (16)is fulfilled.Let us remark that, considering the Hamiltonian at order two in the masses, theBirkhoff normal form is not always convergent at high order, especially when the eccen-tricities are significant or the system is too close to a mean-motion resonance. Indeed,in these cases, the transformation T O , which brings the Hamiltonian at order two inthe masses, induces a big change in the coefficients of the secular Hamiltonian, thatcan prevent the convergence of the normalization procedure. On the contrary, the al-gorithm seems to be convergent at first order in the masses (see the convergence ausens des astronomes in Libert & Henrard (2006)).Assuming that the non-resonance conditions (16) are satisfied up to an order r largeenough, the remainder R ( r ) ( I , ϕ ) can be neglected and we easily obtain an analyticalexpression of the secular frequencies. Indeed, the equations of motion for the truncatedHamiltonian are ˙ I = 0 and ˙ ϕ = ∂H ( r ) ∂ I , (17)and lead immediately to the expression of the two frequencies ˙ ϕ and ˙ ϕ . Let us remarkthat, as the generating functions of the Lie series depend only on the angular difference ϕ − ϕ , the frequency of the apsidal difference ∆̟ = ω − ω is equal to ˙ ϕ − ˙ ϕ . I ( t ) = I (0) and ϕ ( t ) = ϕ (0) + t ˙ ϕ (0) , where I (0) and ϕ (0) correspond to the values of the initial conditions. To validate ourresults, we will compare our analytical integration with the direct numerical integra-tion of the full three-body problem, by using the symplectic integrator SBAB T ( r ) B the canonical transformationinduced by the Birkhoff normalization up to the order r , namely T ( r ) B (cid:0) I , ϕ (cid:1) = exp L X ( r ) ◦ . . . ◦ exp L X (1) (cid:0) I , ϕ (cid:1) . (18)We denote by C ( r ) the composition of all the canonical changes of coordinates definedin Sections 2–4, namely C ( r ) = T F ◦ T O ◦ D ◦ A ◦ T ( r ) B . (19)Taking the initial conditions (cid:0) a (0) , λ (0) , e (0) , ω (0) (cid:1) , we can compute the evolution ofthe orbital elements by using the following scheme (cid:0) a (0) , λ (0) , e (0) , ω (0) (cid:1) (cid:16) C ( r ) (cid:17) − ◦ E − −−−−−→ ( I (0) , ϕ (0)) y(cid:0) a ( t ) , λ ( t ) , e ( t ) , ω ( t ) (cid:1) E ◦ C ( r ) ←−−−−− ( I ( t ) = I (0) , ϕ ( t ) = ϕ (0) + t ˙ ϕ (0)) , (20)where ( Λ , λ , ξ , η ) = E − ( a (0) , λ (0) , e (0) , ω (0)) is the non-canonical change of coor-dinates (2). Thus, the analytical integration via normal form actually reduces to atransformation of the initial conditions to secular action-angles coordinates, the com-putation of the flow at time t in these coordinates, followed by a transformation back tothe original orbital elements. Let us stress that, considering only the secular evolution,the scheme above commutes only if r is equal to infinity.In the following sections, we will compare, for several extrasolar systems, the ana-lytical secular evolution of the eccentricities and apsidal difference with the results ofa direct numerical integration. This kind of comparison has been shown to be a verystressing test (see, e.g., Locatelli & Giorgilli (2007) and Sansottera et al. (2013)) forthe accuracy of the whole algorithm constructing the normal form. E cc en t r i c i t i e s Time (yrs) υ Andromedae numord. 2ord. 1 -60-40-20 0 20 40 600.0e+00 5.0e+03 1.0e+04 1.5e+04 2.0e+04 2.5e+04 3.0e+04 ∆ ϖ ( deg ) Time (yrs)e e numord. 2ord. 1 Fig. 1
Long-term evolution of the eccentricities (top panel) and difference of the longitudesof the pericenters (bottom panel) for the υ Andromedae system ( a /a = 0 . SBAB υ Andromedae system
We aim to investigate the improvements of the secular approximation at order two inthe masses, introduced in the previous sections, in describing the long-term evolutionof extrasolar systems close to a mean-motion resonance. The planetary system υ An-dromedae c and d is well known for his proximity to the 5:1 mean-motion resonance.This has notably been confirmed analytically in Libert & Henrard (2007), where theauthors argued that a first order approximation gives a good qualitative approxima-tion of the secular dynamics of the system. In the following, we show that a secondorder approximation can quantitatively enhance the determination of the secular fre-quencies, as well as the extremal values of the eccentricities and difference in apsidal E cc en t r i c i t i e s Time (yrs) υ Andromedae ( a /a = 0.335 ) numord. 2ord. 1 E cc en t r i c i t i e s Time (yrs) υ Andromedae ( a /a = 0.338 ) e e e e numord. 2ord. 1 Fig. 2
Long-term evolution of slightly different versions of the υ Andromedae system, wherethe semi-major axis of planet d has been modified to be closer to the 5:1 resonance. For a /a = 0 .
335 (left panel), the secular approximation at order two in the masses is stillefficient, while it is no more suitable for a /a = 0 .
338 (right panel). angles reached during the long-term evolution of the planets. For this study, we usethe orbital parameters reported in Wright et al. (2009) † .In order to take into account the proximity of the system to the 5:1 mean-motionresonance, we must include, in the approximation at order two in the masses, the effectsof all the terms up to the trigonometric degree 6 in the fast angles and up to degree4 in the secular variables, namely we set K F = 6 and K S = 4 (see the definitionsin Subsection 3.1). After having constructed the secular approximation, we perform aBirkhoff normal form up to order r = 10 (see Subsection 4.1), which corresponds totaking into account the secular variables up to order 12.In Figure 1, we report the long-term evolution of the eccentricities (top panel)and difference of the longitudes of the pericenters (bottom panel) obtained analyti-cally with our second order approximation (red curves). We compare it to the directnumerical integration of the full three-body problem (i.e., including the fast motions)in heliocentric canonical variables (green curves). The agreement between both curvesis excellent; the second order theory reproduces qualitatively and quantitatively theresults of the numerical integration. A comparison with the first order approximationis also shown (blue curves) and gives evidence of the improvement of the second orderapproximation for systems close to a mean-motion resonance.To highlight the dependency on the truncation parameters, we report, in the tablebelow, the values of the secular period for different values of K F and K S . The periodobtained via numerical integration is ∼ K F = 6 and K S = 4. This validates our choice of the truncation limits. † Let us note that more recent parametrizations consistent with a 30 ◦ mutual inclinationof the two planets (McArthur et al. (2010)) and a fourth planet in the system (Curiel et al.(2010)) have been introduced.3 K F K S Secular period4 2 71326 4 70358 6 6998In Figure 2, we slightly modify the semi-major axis of planet d in such a waythat the modified υ Andromedae systems are closer and closer to the 5:1 mean-motionresonance. First we set a /a = 0 .
335 (left panel). In this case the approximationat order one is not good enough, while the one at order two is still suitable for thecomputation of the long-term evolution of the system, even if the approximation isworst than the one corresponding to the real υ Andromedae in Figure 1. Finally, setting a /a = 0 .
338 (right panel), both the secular approximations completely fail. Indeed,in this case, the system is too close to the resonance to be qualitatively described by asecular approximation.
On the contrary to a first order analytical theory, an expansion to the second order inthe masses takes into account the influence of the initial values of the fast angles onthe secular evolution of the system. Let us stress that, as the averaging process givingthe first order secular Hamiltonian (5) does not involve any canonical transformation,we take as “averaged” initial conditions the original ones † .A change in the mean anomaly of a planet can have a significant impact on thesecular period of the system. To show this, we plot, in Figure 3, the extrasolar sys-tem HD 169830 for two different values of the inner planet mean anomaly: M = 0 ◦ and M = 160 ◦ , all the other orbital parameters being unchanged and issued fromMayor et al. (2004). The displacement between the two secular evolutions is obvious inthe top panel. Let us note that, for both values, our second order averaged Hamiltonian(in red) is very accurate and coincides with the numerical evolution (in green). Thelimitation of the secular expansion to order one in the masses is pointed out in thebottom panel of Figure 3. The first order evolution (blue curve) is the same regardlessthe initial value of the mean anomaly, on the contrary to the approximation at ordertwo in the masses (red curves). We now study the proximity to a mean-motion resonance of the two-planet exosys-tems discovered so far. This represents an extension of the results in Libert & Henrard(2007), previously obtained with an approximation at order one in the masses.Let us make some heuristic considerations. For systems that are very close to a low-order mean-motion resonance k : k , i.e., k n ∗ − k n ∗ ≈
0, the generating functionsrelated to the second order approximation (i.e., µ χ ( O and µ χ ( O defined in (6)and (7), respectively) contain the so-called small divisors . The presence of small divisorsis a major problem in perturbation theory, and here can prevent the convergence of † For sake of completeness, we check that computing the “averaged” initial conditions usingthe generating functions χ and χ , as in the approximation at order two in the masses, doesnot influence neither qualitatively nor quantitatively the results.4 E cc en t r i c i t i e s Time (yrs)
HD 169830 ( M =0 ° vs M =160 ° ) numord. 2 E cc en t r i c i t i e s Time (yrs) ord. 2ord. 1
Fig. 3
Influence of the initial mean anomaly, M , on the secular evolution of the HD 169830system. Top: long-term evolution for M = 0 ◦ and M = 160 ◦ . In both cases, the second orderapproximation (red curves) reproduce accurately the numerical integration (green curves).Bottom: comparison between the evolution at order one (blue curve) and two (red curves) inthe masses. See text for more details. the second order averaging over the fast angles. Instead, for a system that is only near to a mean-motion resonance, but not too close, the approximation at order twoin the masses, including the main effects of the nearest low-order resonance, enablesto describe with great accuracy the long-term evolution of the system. Finally, thesecular evolution of a system that is far from any low order mean-motion resonance isaccurately depicted by the approximation at order one in the masses. Indeed, in thiscase, we can safely replace the canonical transformation T O with the classical firstorder average over the fast angles, see equation (5).Let us go into details. To evaluate the proximity of a planetary system to a mean-motion resonance, we introduce a criterion similar to the one in Libert & Henrard(2007). The idea is to rate the proximity to a mean-motion resonance by looking atthe canonical change of coordinates induced by the approximation at order two inthe masses. Precisely, we consider the low order terms of the canonical transformation induced by T O , writing the averaged variables ( ξ ′ , η ′ ) as ξ ′ j = ξ j − ∂ µ χ ( O ∂η j = ξ j − ξ j ∂ µ χ ( O ∂η j ! ,η ′ j = η j − ∂ µ χ ( O ∂ξ j = η j − η j ∂ µ χ ( O ∂ξ j ! , for j = 1 , µ χ ( O carries the main infor-mation about the proximity to a mean-motion resonance, and we will focus here onthe coefficients of the functions δξ j = 1 ξ j ∂ µ χ ( O ∂η j and δη j = 1 η j ∂ µ χ ( O ∂ξ j . (21)In these expressions, we aim to determine the most important periodic terms whosecorresponding harmonic k · λ identifies the main important mean-motion resonance tothe system. For each system, we define a radius ̺ of a polydisk ∆ ̺ around the originof R , ∆ ̺ = n ( ξ , η ) ∈ R : ξ j + η j ≤ ̺ j , j = 1 , o , so as to include in that domain the initial conditions. Given an analytic function f ,j ( λ , ξ , η ) of the form (4), that reads f ( λ , ξ , η ) = X k f ( k ) ( ξ , η ) sincos ( k · λ ) , where f ( k ) ( ξ , η ) = X | l | + | m | = j f ( k ) l , m ξ l η m , we can easily bound the sup-norm of the terms corresponding to the harmonic k · λ in f , by bounding f ( k ) in the polydisk ∆ ̺ with the norm k f ( k ) k ̺ = X l , m | f ( k ) l , m | ̺ l + m ̺ l + m . (22)Applying the same criterion, for each angular combination k · λ , we evaluate k δξ ( k ) j k ̺ and k δη ( k ) j k ̺ for j = 1 , δξ ∗ j = max k ( k δξ ( k ) j k ̺ ) and δη ∗ j = max k ( k δη ( k ) j k ̺ ) . For convenience we also introduce the following parameters: δ j = min( δξ ∗ j , δη ∗ j ) for j = 1 , δ = max( δ , δ ) . The parameter δ is a measure of the change from theoriginal secular variables to the averaged ones. The actual computation of δ is quitecumbersome, but is more reliable than just looking at the semi-major axes ratio, sinceit holds information about the non-linear character of the system. Table 1
Evaluation of the proximity to a mean-motion resonance (MMR). We report here the values of a /a , µ , k n ∗ + k n ∗ (where k and k correspondto the mean-motion resonance in brackets), δ and δ for each system considered. For our study, we use the following parametrizations: Wright et al.(2009) for HD 11964, HD 12661, υ Andromedae, HD 108874 and HD 183263; Meschiari et al. (2011) for HD 74156 and HD 177830; Jones et al. (2010)for HD 134987; Giguere et al. (2012) for HD 163607; S´egransan et al. (2009) for HD 147018; Tuomi & Kotiranta (2009) for HD 11506; H´ebrard et al.(2010) for HD 9446; Mayor et al. (2004) for HD 169830; JPL at the Julian Date 24404005 for the Sun-Jupiter-Saturn system; Wittenmyer et al. (2009)for HD 128311. See text for more details.System a /a µ k n ∗ + k n ∗ δ δ S ec u l a r HD 11964 0 .
072 5 . × − .
283 (51:1) 5 . × − sin( − λ + λ ) 9 . × − sin( − λ + 2 λ )HD 74156 0 .
075 6 . × − .
579 (48:1) 9 . × − cos( 4 λ − λ ) 3 . × − cos( λ − λ )HD 134987 0 .
140 1 . × − .
052 (19:1) 5 . × − sin( − λ + λ ) 9 . × − sin( − λ + 2 λ )HD 163607 0 .
149 2 . × − .
686 (17:1) 1 . × − cos( 3 λ − λ ) 3 . × − sin( − λ + 2 λ )HD 12661 0 .
287 1 . × − .
671 (6:1) 1 . × − sin( − λ + λ ) 1 . × − sin( − λ + 2 λ )HD 147018 0 .
124 6 . × − .
557 (22:1) 2 . × − sin( − λ + λ ) 1 . × − sin( − λ + 2 λ ) n e a r a MM R HD 11506 0 .
263 2 . × − .
720 (7:1) 2 . × − sin( − λ + 7 λ ) 2 . × − cos( λ − λ )HD 177830 0 .
420 9 . × − .
889 (4:1) 2 . × − cos( λ − λ ) 1 . × − cos( λ − λ )HD 9446 0 .
289 1 . × − .
048 (6:1) 2 . × − sin( − λ + λ ) 2 . × − sin( − λ + 2 λ )HD 169830 0 .
225 2 . × − .
358 (9:1) 1 . × − cos( λ − λ ) 2 . × − cos( λ − λ ) υ Andromedae 0 .
329 3 . × − .
505 (5:1) 1 . × − cos( λ − λ ) 8 . × − cos( λ − λ )Sun-Jup-Sat 0 .
546 9 . × − .
010 (5:2) 1 . × − cos( − λ + 2 λ ) 2 . × − cos(2 λ − λ ) MM R HD 108874 0 .
380 1 . × − .
338 (4:1) 1 . × − cos( λ − λ ) 4 . × − sin( − λ + 4 λ )HD 128311 0 .
622 3 . × − .
924 (2:1) 6 . × − sin( − λ + 2 λ ) 1 . × − sin( − λ + 2 λ )HD 183263 0 .
347 3 . × − .
107 (5:1) 2 . × − cos( λ − λ ) 5 . × − cos( λ − λ )7 HD 108874
HD 128311
HD 183263
HD 169830 υ Andromedae
Sun-Jupiter-Saturn
HD 11506
HD 177830
HD 9446
HD 163607
HD 12661
HD 147018
HD 11964
HD 74156
HD 134987
Fig. 4
Long-term evolution of the eccentricities of the extrasolar systems of Table 1, obtainedby three different ways: (i) the direct numerical integration via
SBAB
The results for all the extrasolar systems we have considered are shown in Table 1and Figure 4. In Table 1, we report, for each system, the numerical values of the semi-major axes ratio a /a , the mass ratio µ , the small divisor k n ∗ + k n ∗ (where k and k correspond to the mean-motion resonance in brackets) and the two aforementionedparameters, δ and δ , that will be used to set up a criterion evaluating the proximityto a mean-motion resonance. In Figure 4, we plot the evolution of the eccentricitiesobtained by direct numerical integration of the Newton equations (green curves) and the ones obtained with a secular Hamiltonian at order one (blue curves) and two (redcurves) in the masses. As in Section 5, we limit the Birkhoff normal form at order r = 10 (i.e., 12 in the secular variables). Let us stress that, for the HD 128311 andHD 183236 systems, due to their close proximity to the mean-motion resonances 2:1and 5:1, respectively, the canonical transformation T O performing the second orderapproximation is not close to the identity. Therefore, we report only their numericalintegrations.Comparing the data in Table 1 and the corresponding plots in Figure 4, we canroughly distinguish three cases: (i) if δ < . × − , the first order approximationdescribes the secular evolution with great accuracy, therefore we label these systemsas secular ; (ii) if 2 . × − < δ ≤ . × − , a second order average of the Hamil-tonian is required to describe the long-term evolution in detail, we label them as nearmean-motion resonance (the υ Andromedae system analyzed in Section 5 is the typ-ical example of such category); (iii) if δ > . × − , the system is too close to amean-motion resonance and a secular approximation is not enough to describe theirevolution, then we label them as in mean-motion resonance . In this case it would beworthwhile to consider a resonant Hamiltonian instead of a secular approximation.Let us note that the Sun-Jupiter-Saturn system and HD 108874 are both really closeto the border between near mean-motion resonance and in mean-motion resonance categories. Indeed, a much refined secular approximation could be used, for instanceincreasing the values of K F and K S , without having to resort to a resonant model.The criterion introduced above is clearly heuristic and quite rough, nevertheless wethink it is useful to discriminate between the different behaviors of planetary systems. In this work we have analyzed the long-term evolution of several exoplanetary systemsby using a secular Hamiltonian at order two in the masses. The second order approx-imation, as explained in detail in Section 3, includes a careful treatment of the maineffects due to the proximity to a low-order mean-motion resonance.Starting from the secular Hamiltonian, we have computed a high-order Birkhoffnormal form via Lie series, introducing action-angle coordinates for the secular vari-ables. This enabled us to compute analytically the evolution on the secular invarianttorus and to obtain the long-term evolution of the eccentricities and apsidal difference.As a result, for all the systems that are not too close to a mean-motion resonance,we have shown an excellent agreement with the direct numerical integration of the fullthree-body problem (including the fast motions). The influence of the mean anomalieson the secular evolution of the systems has also been pointed out. Furthermore, evalu-ating the difference between the original and the averaged secular coordinates, we haveset up a simple (and rough) criterion to discriminate between three different behav-iors: (i) secular system, where a first order approximation is enough; (ii) system near amean-motion resonance, where an approximation at order two is required; (iii) systemthat are really close to or in a mean-motion resonance, where a resonant model shouldbe used. In particular, we find that HD 11964, HD 74156, HD 134987, HD 163607,HD 12661 and HD 147018 belong to (i); HD 11506, HD 177830, HD 9446, HD 169830and υ Andromedae to (ii); HD 108874, HD 128311 and HD 183263 to (iii). Let us remark that these results could be extended to the spatial case with minorchanges. Indeed, after the reduction of the angular momentum, the starting Hamilto-nian would have exactly the same form as H ( T ) , defined in (4).Moreover, having such a good analytical description of the orbits, even for systemsthat are near a mean-motion resonance, we can also study the effective stability ofextrasolar planetary systems in the framework of the KAM and Nekhoroshev theories.This topic deserves further investigation in the future.Finally, a natural extension to the present work would be the study of the secu-lar evolution of systems that are really close to or in a mean-motion resonance. Aspreviously said, a resonant Hamiltonian that keeps the dependency on the resonantcombinations of the fast angles has to be considered. This study is reserved for futurework. Acknowledgements
The work of A.-S. L. is supported by an FNRS Postdoctoral ResearchFellowship. The work of M. S. is supported by an FSR Incoming Post-doctoral Fellowship ofthe Acad´emie universitaire Louvain, co-funded by the Marie Curie Actions of the EuropeanCommission.
A Secular Hamiltonian for the υ Andromedae up to order 6
We report here the expansion of the secular Hamiltonian H (sec) (equation (10)) of the υ Andromedae extrasolar system up to degree 6 in ( ξ , η ) . In particular, we show both theapproximations at order one and two in the masses to highlight the differences. A detaileddescription of the υ Andromedae system is given in Section 5. As this system is near the 5:1mean-motion resonance, the main difference between the two secular approximations affectsterms that are at least of order 6 in the canonical secular variables. ξ ξ η η First order Second order0 0 0 0 − . × +0 − . × +0 − . × − − . × − . × − . × − − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − . × − . × − − . × − − . × − . × − . × − . × − . × − − . × − − . × − . × − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − − . × − . × − . × − . × − − . × − − . × − . × − . × − . × − − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − − . × − . × − . × − − . × − − . × − . × − . × − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − − . × − . × − − . × − − . × − . × − . × − − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − . × − − . × − . × − . × − − . × − − . × − . × − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − . × − − . × − − . × − − . × − . × − . × − − . × − − . × − . × − . × − − . × − − . × − . × − References [1] V.I. Arnold:
Mathematical Methods of Classical Mechanics - Second Edition , Springer, Grad-uate Texts in Mathematics, (1989).[2] C. Beaug´e, D. Nesvor´ny, L. Dones: A high-order analytical model for the secular dynamics ofirregular satellites , AJ, , 2299–2313 (2006).[3] G.D. Birkhoff,
Dynamical systems , AMS Ccolloquium Publications, IX (1927).[4] A. Celletti and L. Chierchia: KAM stability and celestial mechanics , Mem. Amer. Math. Soc., , 1–134 (2007).[5] S. Curiel, J. Canto, L. Georgiev, C. Chavez, A. Poveda:
A fourth planet orbiting upsilonAndromedae , A&A, , A78 (2010).[6] A. Deprit:
Canonical transformations depending on a small parameter , Celestial Mechanics, , 12–30 (1969).[7] J. Fejoz: D´emonstration du “th´eor`eme d’Arnold” sur la stabilit´e du syst`eme plan´etaire(d’apr`es Michael Herman) , Ergodic Theory Dynam. Systems, N.5, 1521–1582 (2005).[8] E.B. Ford, B. Kozinsky, F.A. Rasio:
Secular evolution of hierarchical triple star systems , ApJ, , 385–401 (2000).[9] F. Gabern, A. Jorba, U. Locatelli:
On the construction of the Kolmogorov normal form forthe Trojan asteroids , Nonlinearity, N.4, 1705–1734 (2005).[10] A. Giorgilli:
Quantitative methods in classical perturbation theory , From Newton to chaos:modern techniques for understanding and coping with chaos in N-body dynamical systems,Nato ASI school, A.E. Roy and B.D. Steves eds. (1995).[11] A. Giorgilli and U. Locatelli:
Introduction to canonical perturbation theory for nearly integrablesystems , Chaotic worlds. In: Proceedings of the Nato Advanced Study Institute, (2003).[12] A. Giorgilli, U. Locatelli, M. Sansottera:
Kolmogorov and Nekhoroshev theory for the problemof three bodies , CeMDA, , 159–173 (2009).[13] A. Giorgilli and M, Sansottera:
Methods of algebraic manipulation in perturbation theory ,Workshop Series of the Asociacion Argentina de Astronomia, , 147–183 (2011).[14] M.J. Giguere, D.A. Fischer, A.W. Howard et al.: A high eccentricity component in the doubleplanet system around HD 163607 and a planet around HD 1645091 , ApJ, , id. 4(2012).[15] G. H´ebrard, X. Bonfils, D. S´egransan et al.:
The SOPHIE search for northern extrasolarplanets II. A multi-planet system around HD9446 , A&A, , id. A69 (2010).[16] J. Henrard:
The algorithm of the inverse for Lie transform , Recent Advances in DynamicalAstronomy, Astrophysics and Space Science Library, , 248–257 (1973).[17] G.-I. Hori, Theory of general perturbations with unspecified canonical variables , Publ. Astron.Soc. Jpn., , 287–296 (1966).[18] H.R.A. Jones, R.P. Butler, C.G. Tinney et al.: A long-period planet orbiting a nearby Sun-likestar , MNRAS, , 1703–1713 (2010).[19] B. Katz, S. Dong, R. Malhotra:
Long-Term Cycling of Kozai-Lidov Cycles: Extreme Eccen-tricities and Inclinations Excited by a Distant Eccentric Perturber , PRL, , 181101(2011).[20] K. Kholshevnikov:
D’Alembertian Functions in Celestial Mechanics , Astronomy Reports, The Hamiltonian in the Planetary or Satellite Problem as a d’AlembertianFunction , Astronomy Reports, , 577–579 (2001).[22] A.N. Kolmogorov: Preservation of conditionally periodic movements with small change in theHamilton function , Dokl. Akad. Nauk SSSR, , 527 (1954). Engl. transl. in: Los AlamosScientific Laboratory translation LA-TR-71-67; reprinted in: Lecture Notes in Physics .[23] J. Laskar: Secular evolution over 10 million years , A&A, , 341–362 (1988).[24] J. Laskar:
Syst`emes de variables et ´el´ements , Les M´ethodes modernes de la M´ecanique C´eleste,Editions Fronti`eres, 63–87 (1989).[25] J. Laskar and P. Robutel:
Stability of the Planetary Three-Body Problem — I. Expansion ofthe Planetary Hamiltonian , CeMDA, , 193–217 (1995).[26] J. Laskar and P. Robutel: High order symplectic integrators for perturbed Hamiltonian systems ,CeMDA, , 39–62 (2001).[27] M.H. Lee and S.J. Peale: Secular Evolution of Hierarchical Planetary Systems , ApJ, ,1201–1216 (2003).[28] A.-S. Libert and J. Henrard:
Analytical approach to the secular behaviour of exoplanetarysystems , CeMDA, , 187–200 (2005).2[29] A.-S. Libert and J. Henrard: Secular apsidal configuration of non-resonant exoplanetary sys-tems , Icarus, , 186–192 (2006).[30] A.-S. Libert and J. Henrard:
Analytical study of the proximity of exoplanetary systems tomean-motion resonances , A&A, , 759–763 (2007).[31] A.-S. Libert and N. Delsate:
Interesting dynamics at high mutual inclination in the frameworkof the Kozai problem with an eccentric perturber , MNRAS, , 2725–2736 (2012).[32] U. Locatelli and A. Giorgilli:
Construction of Kolmogorov’s normal form for a planetarysystem , Regular and Chaotic Dynamics , 153–171 (2005) .[33] U. Locatelli and A. Giorgilli: Invariant tori in the Sun-Jupiter-Saturn system , DCDS-B, ,377–398 (2007).[34] M. Mayor, S. Udry, D. Naef et al.: The CORALIE survey for southern extra-solar planetsXII. Orbital solutions for 16 extra-solar planets discovered with CORALIE , A&A, ,391–402 (2004).[35] S. Meschiari, G. Laughlin, S.S. Vogt et al.:
The LICK-CARNEGIE survey: Four new exoplanetcandidates , ApJ, , 117–128 (2011).[36] B.E. McArthur, G.F. Benedict, R. Barnes, E. Martioli, S. Korzennik, Ed. Nelan, R.P. Butler:
New Observational Constraints on the υ Andromedae System with Data from the HubbleSpace Telescope and Hobby-Eberly Telescope , ApJ, , 1203–1220 (2010).[37] S. Naoz, W.M. Farr, Y. Lithwick, F.A. Rasio, J. Teyssandier:
Hot Jupiters from secular planet-planet interactions , Nature, , 187–189 (2011).[38] H. Poincar´e:
Les m´ethodes nouvelles de la M´ecanique C´eleste , Gauthier-Villars (1892).[39] H. Poincar´e:
Le¸cons de M´ecanique C´eleste, tomes I–II , Gauthier-Villars (1905).[40] P. Robutel:
Stability of the Planetary Three-Body Problem — II. KAM Theory and Existenceof Quasiperiodic Motions , CeMDA, , 219–261 (1995).[41] M. Sansottera, U. Locatelli, A. Giorgilli: A Semi-Analytic Algorithm for Constructing LowerDimensional Elliptic Tori in Planetary Systems , CeMDA, , 337–361 (2011).[42] M. Sansottera, U. Locatelli, A. Giorgilli:
On the stability of the secular evolution of the planarSun-Jupiter-Saturn-Uranus system , Math. Comput. Simulat., , 1–14 (2013)[43] D. Segransan, S. Udry, M. Mayor et al.: The CORALIE survey for southern extrasolar planets.XVI. Discovery of a planetary system around HD 147018 and of two long period andmassive planets orbiting HD 171238 and HD 204313 , A&A, , id.A45 (2009).[44] M. Tuomi & S. Kotiranta:
Bayesian analysis of the radial velocities of HD 11506 revealsanother planetary companion , A&A, , L13–L16 (2009).[45] D. Veras and P.J. Armitage:
Extrasolar Planetary Dynamics with a Generalized PlanarLaplace-Lagrange Secular Theory , ApJ, , 1311–1322 (2007).[46] R.A. Wittenmyer, M. Endl, W.D. Cochran et al.:
A search for multi-planet systems using theHOBBY-EBERLY Telescope , ApJSS, , 97–119 (2009).[47] J.T. Wright, S. Upadhyay, G.W. Marcy, D.A. Fisher et al.:
Ten new and updated multiplanetsystems and a survey of exoplanetary systems , ApJ,693