Compact planetary systems perturbed by an inclined companion: I. Vectorial representation of the secular model
CCompact planetary systems perturbed by an inclined companion:I. Vectorial representation of the secular model
Gwenaël Boué , and Daniel C. Fabrycky [email protected] ABSTRACT
The non-resonant secular dynamics of compact planetary systems are modeled by a perturbingfunction which is usually expanded in eccentricity and absolute inclination with respect to theinvariant plane. Here, the expressions are given in a vectorial form which naturally leads to anexpansion in eccentricity and mutual inclination. The two approaches are equivalent in mostcases, but the vectorial one is specially designed for those where a quasi-coplanar system tilts asa whole by a large amount. Moreover, the vectorial expressions of the Hamiltonian and of theequations of motion are slightly simpler than those given in terms of the usual elliptical elements.We also provide the secular perturbing function in vectorial form expanded in semimajor axisratio allowing for arbitrary eccentricities and inclinations. The interaction between the equatorialbulge of a central star and its planets is also provided, as is the relativistic periapse precession ofany planet induced by the central star. We illustrate the use of this representation for followingthe secular oscillations of the terrestrial planets of the solar system, and for Kozai cycles as maytake place in exoplanetary systems.
Subject headings: methods: analytical — methods: numerical — celestial mechanics — planets andsatellites: dynamical evolution and stability — planets and satellites: general — planet-star interactions
1. Introduction
Observations show nearly half of solar-typestars have a stellar companion (Duquennoy &Mayor 1991; Raghavan et al. 2010) and close-ingiant planets are found in such binary systems(Zucker & Mazeh 2002; Udry & Santos 2007).These observations have motivated extensive stud-ies of the dynamics of a single planet evolvingin a binary stellar system (Holman et al. 1997;Wu & Murray 2003; Fabrycky & Tremaine 2007;Wu et al. 2007; Correia et al. 2011; Naoz et al.2011; Lithwick & Naoz 2011; Katz et al. 2011;Veras & Tout 2012; Kratter & Perets 2012; Naozet al. 2012, 2013a,b). But above all, these systems Department of Astronomy and Astrophysics, Univer-sity of Chicago, 5640 South Ellis Avenue, Chicago, IL60637, USA Astronomie et Systèmes Dynamiques, IMCCE-CNRSUMR 8028, Observatoire de Paris, UPMC, 77 Av. Denfert-Rochereau, 75014 Paris, France. present interesting dynamics due to the Lidov-Kozai mechanism (Lidov 1962; Kozai 1962) whichalso provides a natural way to form hot Jupiters.The motion described by Lidov-Kozai dynamicstakes place at large inclination and has the struc-ture of a 1:1 resonance between the precession fre-quency of the longitude of pericenter (cid:36) and of thelongitude of the ascending node Ω of the planet.Inside the resonance, the critical angle, equal tothe argument of the pericenter ω = (cid:36) − Ω , libratesaround either 90 ◦ or 270 ◦ , and the eccentricity andinclination undergo large amplitude oscillations inantiphase. This behavior combined with tidal dis-sipation is able to shrink the orbit of a cold Jupiterdown to orbital periods of about 3 days. Dur-ing the evolution, the apsidal precession becomesdominated by general relativity, in which case thesystem exits the Lidov-Kozai resonance and theorbit evolves through tides as if the companionwas nonexistent. The final obliquity of the starwith respect to the planet orbit is often quite large1 a r X i v : . [ a s t r o - ph . E P ] M a y Fabrycky & Tremaine 2007; Correia et al. 2011;Naoz et al. 2011, 2012; Li et al. 2014; Petrovich2014). These theoretical predictions are supportedby many observations of misaligned systems witha hot Jupiter (Winn et al. 2010; Triaud et al. 2010;Triaud 2011; Albrecht et al. 2012). The dynamicalstructure associated to the Lidov-Kozai mecha-nism subsists when the outer stellar companion isreplaced by a planet (Terquem & Papaloizou 2002;Michtchenko et al. 2006; Libert & Henrard 2007,2008; Migaszewski & Goździewski 2009; Mardling2010; Naoz et al. 2011; Libert & Delsate 2012).The outer planet only needs to be placed on acloser orbit than in the stellar case such that thesecular timescale does not exceed the lifetime ofthe system. This framework has mainly been usedto study the evolution and/or the formation ofJupiters on eccentric orbits. If both planets formin the same protoplanetary disk, the initial mutualinclination is expected to be small, but it can begenerated by high order resonance crossings dur-ing planetary migration (Libert & Tsiganis 2009,2011), secular resonance overlaps (Wu & Lith-wick 2011), or planet-planet scattering (Nagasawa& Ida 2011). Nevertheless, the latter evolutionsare too violent for the Lidov-Kozai mechanism toplay a significant role (Beaugé & Nesvorný 2012).In multiplanet systems surrounded by an outerstellar companion, apsidal precession frequenciesare dictated by the companion and by the planet-planet interactions. As a consequence, even athigh inclination, if the planet system is sufficientlypacked, planet-planet interactions dominate theapsidal motion, the evolution is stabilized withrespect to the Lidov-Kozai mechanism, eccentrici-ties remain small, and all planets move in concert(Innanen et al. 1997; Takeda et al. 2008; Saleh& Rasio 2009). These systems are classified asdynamically rigid. Although the Lidov-Kozai evo-lution is quenched, the planetary mean plane stillprecesses if it is inclined relative to the orbit ofthe companion. The long term evolution of suchsystems can be followed with Laplace-Lagrangesecond order secular theory. This approach hasalready been applied in the context of planet for-mation in a binary system or of satellite formationaround Uranus (Batygin et al. 2011; Batygin 2012;Morbidelli et al. 2012; Lai 2014). But higher ordersecular theories that assume absolute inclinationremains small for all time can no longer follow the system. However, the secular dynamics should notsimply be followed with Laplace-Lagrange secularmodel, which predicts quasiperiodic, bounded ec-centricities and inclinations. For example, thesolar system itself has chaotic secular dynamics(Laskar 1989, 1990), and only theories accurateto fourth or higher orders in eccentricities andinclinations are able to show the appearance ofthis chaos (Laskar 1984; Lithwick & Wu 2011).The purpose of this paper is to provide a set ofequations describing the secular evolution of non-resonant conservative gravitational systems with amassive central body. The formalism is very gen-eral and can be applied to many different typesof systems. In a subsequent paper (Boué andFabrycky 2014), we use it to study analyticallythe evolution of the spin-orbit angle in compactplanetary systems perturbed by an inclined com-panion. Such analyses can a priori be performednumerically using an n -body integrator. How-ever, the huge difference between orbital periodsand secular timescales makes this approach hardlyfeasible in a reasonable amount of time. For in-stance, the 55 Cnc system (Fischer et al. 2008), acompact multiplanet system perturbed by a stel-lar companion (Mugrauer et al. 2006), containsfive planets whose the innermost has a period of0.73 days whereas the precession motion of theplanetary system is about 70 Myr (Kaib et al.2011, fig.1). The integration should thus last sev-eral hundreds of million years with a time step ofa small fraction of a day. The long term evolutioncan also be followed with Gauss’ method wherethe equations of motion, averaged over the meanlongitudes of all planets, are given by double in-tegrals (Touma et al. 2009). In practice, the firstaveraging is analytical while the second has to becomputed numerically. This semi-analytic tech-nique is faster than n -body codes and it can beapplied to a large class of systems, even those withcrossing orbits. Indeed, it does not assume anyconstraints on inclination, eccentricity, nor semi-major axis ratio. However, compact systems withlarge eccentricities or inclinations are likely to beunstable due to resonances overlap. One can thusassume that in most systems which are at leastregular over a few secular timescales, each pairof planets is non-exclusively either hierarchicalor quasi-coplanar with low eccentricities. Thesehypotheses are motivated by statistical studies2f compact exoplanet systems detected by Ke-pler or by radial velocity (e.g., Tremaine & Dong2012; Figueira et al. 2012; Fabrycky et al. 2012;Wu & Lithwick 2013). Within this framework,planet-planet interactions are expanded either insemimajor axis ratio, or in mutual inclination andeccentricity. Hence, this technique is more re-strictive than Gauss’, but both the Hamiltonianand the equations of motion are analytical. Thismethod is thus faster and more appropriate foranalytical studies.For this study, the expressions are given in avectorial form which is independent of any ref-erence frame. This approach, initiated by Mi-lankovitch (1939), has recently been proved veryefficient in the study of cometary motion (Breiter& Ratajczak 2005), of secular spin-orbit evolutionand spin-spin interaction (Boué & Laskar 2006,2009), of the secular three-body problem (Farago& Laskar 2010; Correia et al. 2011), and of the sec-ular evolution of satellites (Tremaine et al. 2009;Tremaine & Yavetz 2013). A detailed historicaldescription of the construction of these variablesand the associated equations of motion is given inRosengren & Scheeres (2014). The Hamiltonianand the equations of motion, derived in Section 2,are expressed analytically by means of expansionsin eccentricity and mutual inclination on the onehand, and in semimajor axis ratio, on the other.The independence of the vectorial equations fromthe reference frame is particularly useful for prob-lems where the plane of a planet system tilts by ahuge angle with respect to a given reference plane.This specific problem is treated in a subsequentpaper (Boué & Fabrycky, 2014). The model istested in Section 3 against numerical integrations.The conclusions are given in the last section.
2. Formalism
Consider a system composed of an arbitrarynumber of massive bodies orbiting a central mass m . For simplicity, the central body is referred toas the central star and the others are called plan-ets. Nevertheless, the formalism is more generaland can be applied to other systems. For exemple,it can model stars orbiting a black hole, or satellitesystems. We utilize Poincaré canonical astrocen-tric variables composed of astrocentric positionsand barycentric velocities (e.g. Laskar & Robu- tel 1995). Ellipses defined by these variables arenot osculating, but for systems with more thanthree bodies, the formalism is simpler than thatinvolving Jacobi coordinates. The canonical as-trocentric variables (positions and conjugate mo-menta) of each body are noted ( r , ˜ r ) , and the or-bital elements ( a, λ, e, (cid:36), I, Ω) represent the semi-major axis, mean longitude, eccentricity, longi-tude of the pericenter, inclination, and longitudeof the ascending node, respectively. All these ele-ments are given with respect to a fixed referenceframe. Whenever we consider two planets, quan-tities associated to the outermost are noted witha prime such as m (cid:48) . We also denote by J themutual inclination of any pairs of planets. In ad-dition to the previous quantities, we also define β = m m/ ( m + m ) , µ = G ( m + m ) where G is theuniversal gravitational constant, and Λ = β √ µa .All quantities are recalled in Tab. 1. This section is devoted to the expansion of theperturbing function in inclination and eccentric-ity. The typical application of such this approx-imation is to model the interaction between twoplanets that are close to each other. In canonicalastrocentric variables, the Hamiltonian describingthe evolution of a compact pair of planets orbitinga central star H close = K close + H I, close + H D, close (1)is the sum of a Keplerian part K close = − µβ a − µ (cid:48) β (cid:48) a (cid:48) , (2)an indirect part H I, close = ˜ r · ˜ r (cid:48) m , (3)and a direct part H D, close = − G mm (cid:48) a (cid:48) a (cid:48) ∆ , (4)where ∆ = (cid:107) r − r (cid:48) (cid:107) is the distance between thetwo planets. Because we only focus on non-resonant secular evolutions, the indirect part H I, close of the Hamiltonian cancels and the Ke-plerian part K close remains constant. The secularevolution is thus solely controlled by the directpart H D, close .3able 1: Notation.variable Ref. description p l a n e t a r y o r b i t a l a ndph y s i c a l e l e m e n t s r , r (cid:48) astrocentric position ˜ r , ˜ r (cid:48) barycentric velocity a , a (cid:48) semimajor axis α semimajor axis ratio a/a (cid:48) λ , λ (cid:48) mean longitude e , e (cid:48) eccentricity ω , ω (cid:48) argument of pericenter (cid:36) , (cid:36) (cid:48) longitude of pericenter I , I (cid:48) absolute inclination J mutual inclination between two planets ρ sin( J/ σ cos( J/ , Ω (cid:48) longitude of the ascending node ON τ ON + NGO
Fig. 1 origin of longitude N , N (cid:48) Fig. 1 ascending node relative to the reference plane G , G (cid:48) Fig. 1 ascending node relative to another orbit plane ∆ mutual distance (cid:107) r − r (cid:48) (cid:107) m , m (cid:48) planet mass β , β (cid:48) reduced mass µ , µ (cid:48) G ( m + m )Λ , Λ (cid:48) β √ µa e , e (cid:48) eccentricity vector j , j (cid:48) dimensionless angular momentum √ − e ww , w (cid:48) unit vector along the orbital angular momentum ξ , ξ (cid:48) Eq. (15) Souriau variable j + eη , η (cid:48) Eq. (15) Souriau variable j − e T j , V j , W j Eq. (18) Abdullah variables quadratic in inclination and eccentricity s t e ll a r p a r a m e t e r s m star mass R star radius C moment of inertia along the short axis k second fluid Love number J Eq. (36) quadrupole gravitational harmonic ω rotation vector s spin axis L angular momentum Cω s H a m il t o n i a n s H close Eq. (1) Hamiltonian of packed system ¯ H close Eq. (26) secular Hamiltonian of packed system ¯ H hierar Eq. (33) secular Hamiltonian of hierarchical system ¯ H spin Eq. (34) secular Hamiltonian of spin-orbit interaction ¯ H relat Eq. (37) secular Hamiltonian of general relativity b ( k ) s ( α ) Eq. (8) Laplace coefficient f j Eq. (5) coefficients of Le Verrier’s expansion of the perturbing function g j Eq. (25) coefficients of the perturbing function in Abdullah variables c j Tab. 3 coefficients of the perturbing function in Milankovitch variables G gravitational constant c speed of light 4 .1.1. Le Verrier’s expansion The secular component of the expansion of a (cid:48) / ∆ in eccentricity and absolute inclination iswell known and is expressed analytically in, e.g.,(Laskar & Robutel 1995; Ellis & Murray 2000).Although this approach is very convenient becausethe resulting expression is a polynomial in thecanonical Poincaré variables (Laskar & Robutel1995), it has not been designed to study quasi-coplanar systems with large absolute inclinationssuch as compact planetary systems perturbed by adistant and inclined stellar companion. One needsinstead an expansion in mutual inclination as de-rived by Le Verrier (1855). The latter is simplerand more compact than expansions in absolute in-clination because it is a special case where one ofthe absolute inclinations is set to zero. However,the equations of motion given in mutual inclina-tion are more cumbersome, especially in systemswith more than two planets. In the following, werecall the development of the secular componentof a (cid:48) / ∆ exact up to the fourth order in eccentric-ity e, e (cid:48) and mutual inclination J . Then, we re-call Le Verrier’s equations of motion. The secularterms of a (cid:48) / ∆ are (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) = f + f ( e + e (cid:48) − ρ )+ f ee (cid:48) cos( ω − ω (cid:48) )+ f ( e e (cid:48) − ρ ( e + e (cid:48) ))+ f ρ ee (cid:48) cos( ω + ω (cid:48) )+ f ( e + 8 ρ e (cid:48) cos 2 ω (cid:48) )+ f ( e (cid:48) + 8 ρ e cos 2 ω )+ f ρ ee (cid:48) cos( ω − ω (cid:48) )+ f e e (cid:48) cos( ω − ω (cid:48) )+ f ee (cid:48) cos( ω − ω (cid:48) )+ f e e (cid:48) cos 2( ω − ω (cid:48) )+ f ρ , (5)with ρ = sin( J/ . The secular Hamiltonian de-scribing the evolution of a two planet system isthen ¯ H close = − G mm (cid:48) a (cid:48) (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) (6) G’ yx z GO NN’ outer orbitinner orbitreference plane
Fig. 1.— Orbit orientation. ON’N defines thereference plane. O is the origin of longitudes, Nand N’ are the ascending nodes of the orbits of m and m (cid:48) relative to the reference plane, respectively.G and G’ are the ascending nodes of the orbit m and m (cid:48) relative to the orbit m (cid:48) and m , respectively.In (5), the orientation of the orbits are given rel-ative to each other. More precisely, let G be theascending node of the inner orbit relative to theouter one, and similarly, G’ the ascending nodeof the outer orbit with respect to the inner one.G and G’ are thus on the intersection betweenthe two orbit planes, but in opposite directions(see Fig. 1). The angles ω and ω (cid:48) are the ar-guments of pericenter of the two planets relativeto G and G’, respectively. With the notationof Fig. 1, the longitudes of the ascending nodesof m and m (cid:48) relative to the reference plane are Ω = ON and Ω (cid:48) = ON (cid:48) , respectively. FollowingLe Verrier (1855), we denote τ = ON + NG and τ (cid:48) = ON (cid:48) +N (cid:48) G (cid:48) . The two arguments of pericenter ω and ω (cid:48) are then given by ω = (cid:36) − τω (cid:48) = (cid:36) (cid:48) − τ (cid:48) . (7)Functions ( f k ) k =1 ,..., contain the dependency inthe semimajor axis ratio α = a/a (cid:48) . Le Verrier(1855) computed them in terms of Laplace coeffi-cients b ( k ) s ( α ) , which are given by (e.g., Laskar &Robutel 1995)
12 b ( k ) s ( α ) = ( s ) k k ! α k F ( s, s + k, k + 1; α ) , (8)and of their derivatives. In (8), ( s ) = 1 and ( s ) k = s ( s + 1) · · · ( s + k − if k > . The func-tion F ( a, b, c ; x ) is the Gauss hypergeometric func-tion. Le Verrier’s expressions of the functions f k f = 12 b (0)1 / ,f = 18 α b (1)3 / ,f = − α b (0)3 / + 12 (1 + α )b (1)3 / ,f = 932 α b (0)5 / ,f = 98 α b (1)5 / ,f = − α b (0)5 / + 364 α (1 + 3 α )b (1)5 / ,f = − α b (0)5 / + 364 α (3 + α )b (1)5 / ,f = − α (1 + α )b (0)5 / + 34 (2 + 3 α + 2 α )b (1)5 / ,f = − α b (0)5 / + 332 (4 + 9 α )b (1)5 / ,f = − α b (0)5 / + 332 α (9 + 4 α )b (1)5 / ,f = 4564 α b (0)5 / − α (1 + α )b (1)5 / ,f = 218 α b (0)5 / − α (1 + α )b (1)5 / . (9)In the averaged problem, the mean longitudes ( λ, λ (cid:48) ) do not appear anymore, we thus discardtheir evolution. Furthermore, the semimajor axesare constant. The equations of motion derived byLe Verrier (1855) of the other orbital elements ofany two planets m and m (cid:48) whatever are their po- sition relative to each other, are dedt = √ − e Λ e ∂ ¯ H close ∂ω ,d(cid:36)dt = − √ − e Λ e ∂ ¯ H close ∂e + tan I I d Ω dt ,dIdt = sin( τ − Ω)Λ √ − e ∂ ¯ H close ∂J − cos( τ − Ω)Λ √ − e (cid:32) f τ − ρ (cid:112) − ρ ∂ ¯ H close ∂ω (cid:33) ,d Ω dt = − I (cid:32) cos( τ − Ω)Λ √ − e ∂ ¯ H close ∂J + sin( τ − Ω)Λ √ − e (cid:32) f τ − ρ (cid:112) − ρ ∂ ¯ H close ∂ω (cid:33) (cid:33) (10)with f τ = 1sin J (cid:18) ∂ ¯ H close ∂ω + ∂ ¯ H close ∂ω (cid:48) (cid:19) . (11)The relations between ( J, τ ) and the orbital ele-ments relative to the reference plane are cos J = cos I cos I (cid:48) + sin I sin I (cid:48) cos(Ω − Ω (cid:48) ) , sin( τ − Ω) = sin I (cid:48) sin J sin(Ω − Ω (cid:48) ) , cos( τ − Ω) = cos I (cid:48) − cos I cos J sin I sin J . (12)The above equations are sufficient to compute nu-merically the secular evolution of a p -planet sys-tem with low mutual inclinations and low eccen-tricities. However, they have singularities whenany inclination or eccentricity is nil. This is-sue can be circumvented by choosing non singu-lar variables such as k = e cos (cid:36) , h = e sin (cid:36) , q = sin( I/
2) cos Ω , and p = sin( I/
2) sin Ω , andby imposing ˙ p = ˙ q = 0 whenever sin J = 0 . In anycase, the equations of motion (10) are complicatedand do not facilitate the comprehension of the dy-namical behavior of the system. For that reason,we shall consider a new set of variables. The orientation in space and the shape of a Ke-plerian orbit (an ellipse) can be parametrized by6wo vectors: the dimensionless angular momen-tum j = √ − e w and the Laplace-Runge-Lenz(or simply eccentricity) vector e = e u , where w isthe unit vector normal to the orbit along the angu-lar momentum, and u is the unit vector pointingtoward the pericenter. These vectors are used inplace of the usual elliptical elements ( e, (cid:36), I, Ω) .Since each vector has three coordinates, the num-ber of variables increases from four to six. This im-plies that the new variables are not independent,and indeed, they are related by two equations: e · j = 0 , (cid:107) e (cid:107) + (cid:107) j (cid:107) = 1 . (13)The equations of motion expressed in terms ofthese two vectors were derived by Milankovitch(1939). In the secular problem, they read as (e.g.,Breiter & Ratajczak 2005; Tremaine et al. 2009;Rosengren & Scheeres 2014) d j dt = −
1Λ ( j × ∇ j ¯ H + e × ∇ e ¯ H ) ,d e dt = −
1Λ ( e × ∇ j ¯ H + j × ∇ e ¯ H ) , (14)where ¯ H (here ¯ H = ¯ H close ) is the secular Hamilto-nian of the system written in terms of the vectors( e , j , e (cid:48) , j (cid:48) ). The variables ( e , j ) are not singular.Furthermore, they lead to more compact and sym-metrical equations to describe the evolution of thesystem. The expression of the Hamiltonian as afunction of these vectors is presented in the fol-lowing section. This section reproduces the computation of theperturbing function (6) made by Abdullah (2001)in terms of the vectors e and j . This methodinvolves new variables named after Souriau (1969).These variables noted ξ and η , are defined by ξ = j + e , η = j − e . (15)It can easily be shown that (cid:107) ξ (cid:107) = 1 (cid:107) η (cid:107) = 1 . (16) Our expressions differ slightly from those obtained by Ab-dullah (2001) because we use ρ = sin( J/ like in thework of Le Verrier (1855), whereas in Abdullah’s notation, ρ = sin J . Furthermore, the scaling factor in Laplace co-efficients is arbitrary, and Abdullah (2001) do not includethe factor / that we have in (8). Relations (16) let Souriau (1969) conclude that theset of ellipses with one fixed focus and semimajoraxis is equivalent to the product of two spheres: S × S . The dimension of the problem is thusclearly four. From (14) and (15), the derivation ofthe equations of motion of ξ and η is straightfor-ward. The result is d ξ dt = − ξ × ∇ ξ H ,d η dt = − η × ∇ η H . (17)In these variables, the equations of motion (17)are very simple and symmetric. To get the ex-pression of the secular Hamiltonian as a functionof Souriau’s variables, one needs to expand (cid:104) a (cid:48) / ∆ (cid:105) in terms of these variables, but since their normsare equal to one, by construction, they are notsmall quantities. To solve that issue, Abdullah(2001) considered a new family of variables basedon Souriau’s one, and defined by T = 12 (1 − ξ · η ) ,T = 12 (1 − ξ (cid:48) · η (cid:48) ) ,V = 12 (1 − ξ · ξ (cid:48) ) ,V = 12 (1 − η · η (cid:48) ) ,W = 12 (1 − ξ · η (cid:48) ) ,W = 12 (1 − ξ (cid:48) · η ) . (18)These variables can be interpreted as the square ofthe semi-distances between the points representedon the unit sphere by the vectors ξ , ξ (cid:48) , η , and η (cid:48) .Indeed, for instance, T = (cid:18) (cid:107) ξ − η (cid:107) (cid:19) . (19)The variables (18) are quadratic in eccentricitiesand in mutual inclinations. The method is thento express all the quantities appearing in the ex-pression (6), i.e., ρ , e , e (cid:48) , p ( k, (cid:96) ) and q ( k, (cid:96) ) asa function of the variables (18), where p ( k, (cid:96) ) and q ( k, (cid:96) ) are defined for all ( k, (cid:96) ) ∈ N by p ( k, (cid:96) ) = ρ k + (cid:96) e k e (cid:48) (cid:96) cos( kω + (cid:96)ω (cid:48) ) ,q ( k, (cid:96) ) = ρ | k − (cid:96) | e k e (cid:48) (cid:96) cos( kω − (cid:96)ω (cid:48) ) , (20)7f k + (cid:96) is even, and p ( k, (cid:96) ) = ρ k + (cid:96) e k e (cid:48) (cid:96) sin( kω + (cid:96)ω (cid:48) ) ,q ( k, (cid:96) ) = ρ | k − (cid:96) | e k e (cid:48) (cid:96) sin( kω − (cid:96)ω (cid:48) ) , (21)if k + (cid:96) is odd. To simplify the substitution, Abdul-lah (2001) provides recurrence relations to com-pute the p ( k, (cid:96) ) and q ( k, (cid:96) ) which are equivalentto those of Table 2. To initiate the recurrence, wehave e = T ,e (cid:48) = T ,ρ = 12 − − V − V − W − W (cid:112) (1 − T )(1 − T ) ,p (1 ,
0) = V − V + W − W (cid:112) (1 − T )(1 − ρ ) ,p (0 ,
1) = V − V + W − W (cid:112) (1 − T )(1 − ρ ) ,q (1 ,
1) = V + V − W − W p (1 , p (0 , ,p (1 ,
1) = ρ V + V − W − W − − ρ ) p (1 , p (0 , , (22)and q ( k,
0) = p ( k, ,q (0 , (cid:96) ) = − p (0 , (cid:96) ) (23)for all ( k, (cid:96) ) in N . After substituting the terms p ( k, (cid:96) ) and q ( k, (cid:96) ) in (5) by T , T , V , V , W , and W , and truncating at the second order in these new variables, Abdullah (2001) obtained (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) = f + f (2 T + 2 T − V − V − W − W )+ 12 f ( V + V − W − W )+ g ( α T + T )+ g ( V + V )+ g ( W + W )+ g ( T T + V V + W W )+ g ( αT − T )( V + V )+ g ( αT + T )( W + W )+ g ( V W + V W )+ g ( V W + V W ) , (24)where f , f , and f are given in (9) and g = 932 α b (1)5 / ,g = 316 α ( − α − α )b (0)5 / + 332 (2 − α + 2 α ) b (1)5 / ,g = 316 α (5 + 4 α + 5 α )b (0)5 / −
332 (2 + α + 2 α ) b (1)5 / ,g = 916 α b (0)5 / ,g = − α (1 − α )b (0)5 / + 316 (1 − α )b (1)5 / ,g = − α (1 + α )b (0)5 / + 316 (1 + α )b (1)5 / ,g = − α (cid:16) b (0)5 / − α b (1)5 / (cid:17) ,g = − α (cid:16) α b (0)5 / − b (1)5 / (cid:17) . (25)The expansion (24) can be seen as a function ofSouriau’s variables ( ξ , η , ξ (cid:48) , η (cid:48) ) using (18), or asa function of Milankovitch’s variables ( e , j , e (cid:48) , j (cid:48) ) using the definition of Souriau’s variables (15).8able 2: Recurrence relations for the computation of p ( k, (cid:96) ) and q ( k, (cid:96) ) , Eqs. (20)-(21). p ( k, (cid:96) ) = 2 p ( k − , (cid:96) − p (1 , − ρ e e (cid:48) p ( k − , (cid:96) − k ≥ and (cid:96) ≥ p ( k,
1) = 2 p ( k − , p (1 , − ρ e q ( k − , k > p (1 , (cid:96) ) = 2 p (0 , (cid:96) − p (1 ,
1) + ( − (cid:96) ρ e (cid:48) q (1 , (cid:96) − (cid:96) > p ( k,
0) = 2( − k +1 p ( k − , p (1 ,
0) + ρ e p ( k − , k ≥ p (0 , (cid:96) ) = 2( − (cid:96) +1 p (0 , (cid:96) − p (0 ,
1) + ρ e (cid:48) p (0 , (cid:96) − (cid:96) ≥ p (2 ,
1) = 2 p (1 , p (1 , − ρ e q (0 , p (1 ,
2) = 2 p (0 , p (1 ,
1) + ρ e (cid:48) q (1 , q ( k, (cid:96) ) = 2 q ( k − , (cid:96) − q (1 , − e e (cid:48) q ( k − , (cid:96) − k ≥ and (cid:96) ≥ q ( k,
1) = 2 q ( k − , q (1 , − e p ( k − , k ≥ q (1 , (cid:96) ) = 2 q (0 , (cid:96) − q (1 ,
1) + ( − (cid:96) e (cid:48) p (1 , (cid:96) − (cid:96) ≥ q ( k,
0) = p ( k, , q (0 , (cid:96) ) = ( − (cid:96) p (0 , (cid:96) ) ∀ k, l One gets (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) = c + c ( e + e (cid:48) + j · j (cid:48) − c ( e · e (cid:48) )+ c e e (cid:48) + c ( α e + e (cid:48) )+ c ( e · e (cid:48) ) + c (1 − j · j (cid:48) ) + c ( e · j (cid:48) ) + c ( j · e (cid:48) ) + c (cid:0) α (1 − j · j (cid:48) ) e + ( e · e (cid:48) ) e (cid:48) (cid:1) + c (cid:0) (1 − j · j (cid:48) ) e (cid:48) + α ( e · e (cid:48) ) e (cid:1) + c (cid:0) (1 − j · j (cid:48) )( e · e (cid:48) ) − ( e · j (cid:48) )( j · e (cid:48) ) (cid:1) . (26)The first three lines with coefficients c , c , and c correspond to the second order expansion in ec-centricity and mutual inclination. All the otherterms are associated to the order 4. Because thedot product of two vectors is invariant by any rota-tion applied on both vectors, the expression (26)is clearly invariant by rotation. Thus, the refer- ence frame does not need to be aligned with themean plane of the planet system. The coefficients ( c k ) k =1 , ··· , are c = f ,c = 2 f ,c = − f ,c = g ,c = g ,c = 12 ( g + g + g − g − g ) ,c = 12 ( g + g + g + g + g ) ,c = 12 ( g + g − g + g − g ) ,c = 12 ( g + g − g − g + g ) ,c = g + g ,c = g − g ,c = g − g . (27)Their explicit expressions in terms of Laplace coef-ficients are given in Tab.3. The equations of mo-tion, derived from (6), (14), and (26), are writtenin Appendix A.1. We stress here that the vecto-rial expansion (26) is a generalization of the morestandard development in eccentricity and absolute9able 3: Coefficients of the secular expansion ofthe perturbing function in eccentricity and mutualinclination. c = 12 b (0)1 / c = 14 α b (1)3 / c = 34 α b (0)3 / −
12 (1 + α )b (1)3 / c = 916 α b (0)5 / c = 932 α b (1)5 / c = 4532 α b (0)5 / − α (1 + α )b (1)5 / c = 2132 α b (0)5 / − α (1 + α )b (1)5 / c = 1532 α b (0)5 / − α (3 + α )b (1)5 / c = 1532 α b (0)5 / − α (1 + 3 α )b (1)5 / c = − α b (0)5 / + 38 b (1)5 / c = − α b (0)5 / + 38 α b (1)5 / c = 158 α (1 + α )b (0)5 / −
316 (4 + 9 α + 4 α )b (1)5 / inclination (e.g., Laskar & Robutel 1995; Ellis &Murray 2000). Indeed, at low inclination with re-spect to the reference frame, both formalisms areequivalent. If the planet system is tilted by alarge angle, the linear equations of the Laplace-Lagrange approximation are still valid. But athigher orders, where inclinations get coupled witheccentricities, tilted systems can only be describedwith expansions in mutual inclination such as inthe vectorial approach. Consider the case where the two planets m , m (cid:48) are very distant from each other ( α ≡ a/a (cid:48) (cid:46) . ).The Hamiltonian ¯ H hierar governing the secularevolution is similar to ¯ H close (6), ¯ H hierar = − G mm (cid:48) a (cid:48) (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) . (28)The only difference is that (cid:104) a (cid:48) / ∆ (cid:105) is expanded insemimajor axis ratio α rather than in eccentric-ity and inclination. The most important contribu-tions to this development have been made in theXIXth century (Hansen 1853; Hill 1875; Tisserand1889). The explicit expression of the secular per-turbing function expanded at the octupolar order10s (e.g., Laskar & Boué 2010) (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) = 1 + (cid:32) (cid:18) − ρ + 32 ρ (cid:19) (cid:18) e (cid:19) + 154 ρ σ e cos(2 ω ) (cid:33) α (1 − e (cid:48) ) / + 1516 (cid:32) (cid:18) e (cid:19) × (cid:104) σ (cid:0) − ρ + 15 ρ (cid:1) cos( ω − ω (cid:48) )+ ρ (cid:0) − ρ + 15 ρ (cid:1) cos( ω + ω (cid:48) ) (cid:105) + 354 e ρ σ (cid:104) ρ cos(3 ω + ω (cid:48) )+ σ cos(3 ω − ω (cid:48) ) (cid:105)(cid:33) ee (cid:48) α (1 − e (cid:48) ) / , (29)with σ = 1 − ρ . As in section 2.1, one can usethe equations of motion (10) derived by Le Verrier(1855) to get the evolution of the system describedby (29), but once again, the expressions are muchsimpler in a vectorial form. To get the vectorial ex-pression of (cid:104) a (cid:48) / ∆ (cid:105) expanded in semimajor axis ra-tio, we follow an algorithm very similar to the oneof the section 2.1. We define P ( k, (cid:96) ) and Q ( k, (cid:96) ) , ( k, (cid:96) ) ∈ N , as P ( k, (cid:96) ) = σ | k − (cid:96) | p ( k, (cid:96) ) ,Q ( k, (cid:96) ) = σ k + (cid:96) q ( k, (cid:96) ) , (30)where p ( k, (cid:96) ) and q ( k, (cid:96) ) are given in Eqs (20)-(21).The secular part of the expansion in semimajoraxis ratio of the perturbing function is a polyno-mial in α , ρ , e , e (cid:48) , (1 − e (cid:48) ) − / , P ( k, (cid:96) ) , and Q ( k, (cid:96) ) (Tisserand 1889). The recurrence relationssatisfied by P ( k, (cid:96) ) and Q ( k, (cid:96) ) are displayed in Ta-ble 4. The initialization of the recurrence can be In this paper, ω (cid:48) is defined with respect to the ascendingnode G’ of the orbit m (cid:48) relative to the orbit m , while inLaskar & Boué (2010), ω (cid:48) is defined with respect to theascending node G of the orbit m relative to the orbit m (cid:48) .Thus, the two arguments of periastron differ by π . done either in terms of Souriau’s variables e = T ,e (cid:48) = T ,ρ = 12 − − V − V − W − W (cid:112) (1 − T )(1 − T ) ,P (1 ,
0) = V − V + W − W √ − T ,P (0 ,
1) = V − V + W − W √ − T ,Q (1 ,
1) = σ V + V − W − W P (1 , P (0 , ,P (1 ,
1) = ρ V + V − W − W − P (1 , P (0 , , (31)or directly in terms of Milankovitch’s vectors ρ = 1 − w · w (cid:48) , σ = 1 + w · w (cid:48) , P (1 ,
0) = ( w (cid:48) · e ) , P (0 ,
1) = ( w · e (cid:48) ) , Q (1 ,
1) = − σ ( e · e (cid:48) ) + ( e · w (cid:48) )( w · e (cid:48) ) , P (1 ,
1) = − ρ ( e · e (cid:48) ) − ( e · w (cid:48) )( w · e (cid:48) ) , (32)where w = j / √ − e , and w (cid:48) = j (cid:48) / √ − e (cid:48) . Aquick comparison of (31) and (32) suggests thatMilankovitch’s formalism is more adapted for thisproblem. In these variables (cid:104) a (cid:48) / ∆ (cid:105) , expanded insemimajor ratio, reads (cid:28) a (cid:48) ∆ (cid:29) λ,λ (cid:48) = 1 + α j (cid:48) (cid:16) j · j (cid:48) ) − (1 − e ) j (cid:48) − e · j (cid:48) ) (cid:17) + 15 α j (cid:48) (cid:16) ( e · e (cid:48) ) × (cid:2) (1 − e ) j (cid:48) + 35( e · j (cid:48) ) − j · j (cid:48) ) (cid:3) − e · j (cid:48) )( j · e (cid:48) )( j · j (cid:48) ) (cid:17) . (33)One can easily check that (33) is invariant by ro-tation. The associated equations of motion of theplanet and the companion are deduced from (14),(28), and (33) (see Appendix A.2).11able 4: Recurrence relations for the computation of P ( k, (cid:96) ) and Q ( k, (cid:96) ) , Eq. (30). P ( k, (cid:96) ) = 2 P ( k − , (cid:96) − P (1 , − ρ e e (cid:48) P ( k − , (cid:96) − k ≥ and (cid:96) ≥ P ( k,
1) = 2 P ( k − , P (1 , − ρ e Q ( k − , k > P (1 , (cid:96) ) = 2 P (0 , (cid:96) − P (1 ,
1) + ( − (cid:96) ρ e (cid:48) Q (1 , (cid:96) − (cid:96) > P ( k,
0) = 2( − k +1 P ( k − , P (1 ,
0) + σ ρ e P ( k − , k ≥ P (0 , (cid:96) ) = 2( − (cid:96) +1 P (0 , (cid:96) − P (0 ,
1) + σ ρ e (cid:48) P (0 , (cid:96) − (cid:96) ≥ P (2 ,
1) = 2 P (1 , P (1 , − ρ e Q (0 , P (1 ,
2) = 2 P (0 , P (1 ,
1) + ρ e (cid:48) Q (1 , Q ( k, (cid:96) ) = 2 Q ( k − , (cid:96) − Q (1 , − σ e e (cid:48) Q ( k − , (cid:96) − k > and (cid:96) ≥ Q ( k,
1) = 2 Q ( k − , Q (1 , − σ e P ( k − , k > Q (1 , (cid:96) ) = 2 Q (0 , (cid:96) − Q (1 ,
1) + ( − (cid:96) σ e (cid:48) P (1 , (cid:96) − (cid:96) > Q ( k,
0) = P ( k, , Q (0 , (cid:96) ) = ( − (cid:96) P (0 , (cid:96) ) ∀ k, lQ (2 ,
1) = 2 Q (1 , Q (1 , − σ e P (0 , Q (1 ,
2) = 2 Q (0 , Q (1 ,
1) + σ e (cid:48) P (1 , Because of their proper rotations, stars are notspherical and exert a torque on the orbital motionof their planets. Let a system composed of a star m with one planet m . We note s the unit vectoralong the spin axis of the star, J its quadrupolegravitational harmonic, and R its equatorial ra-dius. We make the assumption that s also corre-sponds to the stellar axis of maximal inertia (gy-roscopic approximation). The secular quadrupolepotential energy due to the oblateness of the staracting on the planet is (e.g., BL06) ¯ H spin = G m mJ R a (1 − e ) / (cid:0) − s · w ) (cid:1) . (34)This expression is valid as long as the distance ofthe planet to the star is much larger than the stel-lar radius ( r (cid:29) R ). Nevertheless, since stellardeformations are usually small, we assume that(34) is valid even for close-in planets with verysmall semimajor axis. There is no assumption re-garding the obliquity of the star relative to theorbital plane of the planet. In a generic problemwhere ¯ H represents the Hamiltonian of the sys-tem, the equation of motion satisfied by the spin axis is (BL06) d s dt = − L s × ∇ s ¯ H , (35)where L = Cω s is the angular momentum of thestar, C its moment of inertia along the s axis, and ω its rotation rate. The explicit equations of mo-tion of s , j and e are given in Appendix A.3. Asshown in BL06, the kinetic energy associated tothe rotation of a rigid body around its spin axis,which should be added in the Hamiltonian, doesnot contribute to the equations of motion once theHamiltonian is averaged over the proper rotationof the rigid body, or when the body has an axialsymmetry. As a consequence, we drop this kineticenergy from our equations. The gravity field co-efficient J is deduced from the rotation speed ofthe star according to (e.g., Lambeck 1988) J = k ω R G m , (36)where k is the second Love number. The last effect taken into account is the secu-lar contribution of general relativity on the preces-sion motion of the planets pericenter. For a planet12 orbiting a star m , the associated Hamiltonianreads (e.g., Touma et al. 2009) ¯ H relat = − µ βa c √ − e , (37)where c is the speed of light, β = m m/ ( m + m ) ,and µ = G ( m + m ) . The associated equations ofmotion are given in Appendix A.4. The formalism described above can be used tomodel many different types of systems. For exam-ple, consider a compact planet system perturbedby an outer stellar companion such as the 55 Can-cri system (Kaib et al. 2011). If the planet sys-tem is similar to those detected by the
Kepler spacecraft, it should be dynamically cold with loweccentricities and mutual inclinations. We fur-ther assume that the dynamics of the system isnot dominated by mean-motion resonances. Con-versely, the binary component can be highly ec-centric with large inclination with respect to theplanet plane. In this case, the Hamiltonian of thesystem can be approximated by H tot = (cid:88) ≤ j 3. Numerical tests and applications To check our secular models, we performed sev-eral tests. First, we considered a system composedof the four inner planets of our solar system. We made this choice because the eccentricities and themutual inclinations are well known, and relativelylow. We used the initial conditions provided inYoder (1995). This system is used to test the ex-pansions (6) and (26) of the perturbing functionin small eccentricity and mutual inclination. In asecond step, we chose a hierarchical system under-going Lidov-Kozai oscillations (Lidov 1962; Kozai1962) to test our expansions in semimajor axis (29)and (33). In both cases, the secular evolutions arecompared to numerical simulations done with afull n -body symplectic integrator. Figure 2 displays the evolution over one mil-lion years of the eccentricities and inclinations ofthe four planets of our first system with respectto the initial ecliptic plane ( I Earth = 0 at t = 0 )obtained with an n -body integrator and by solv-ing the Hamiltonian (26). General relativity is in-cluded in both simulations as is the effect of thesolar oblateness. The two integrations give verysimilar results and their distinction is hardly per-ceptible in most subfigures. It is best seen in theevolution of the eccentricity of Venus and of theEarth, where it consists mostly in a slight shift inthe precession frequencies. Thus, the secular for-mulation preserves the main dynamical featuresof this system, as already observed by, e.g., Laskar(2008). The integration of the Hamiltonian (6) us-ing Le Verrier’s equations of motion is not shownhere, because it cannot be distinguished from thesolution obtained with the vectorial approach.In Fig. 3 and 4, we test the effect of a rotation ofthe whole system by an angle ∆ I ∈ { 0, 10, 60, 90,135, 180 } degrees. More precisely, we apply eachof these rotations on the secular integrations only,and we keep the n -body integration of the previ-ous figure unchanged. Since the choice of the ref-erence frame is arbitrary, the eccentricities shouldnot be affected by these rotations. Figure 3 showsthat this is indeed the case whether the evolutionis obtained using Le Verrier’s equations (10) orwith the vectorial approach (14). On the otherhand, Hamiltonians expanded in absolute inclina-tions are not designed to describe the evolution ofsuch highly inclined systems. It is thus natural toobserve discrepancies above ∆ I = 90 ◦ between thesecular model taken from, e.g., Laskar & Robutel(1995) (dotted curve in Fig. 3) and the n -body13 ars time (Myr) time (Myr) inclination ( ◦ ) n -bodyMercury eccentricity Fig. 2.— Comparison between an n -body inte-gration and a secular integration on a system com-posed of the four inner planets of our solar systemfollowed over one million years.integration (solid curve). A similar result is ob-served on Mercury’s inclination in Fig. 4. To makethis figure, we integrated the system in a tiltedframe, and then we applied a rotation of − ∆ I onthe output to place the system back in the refer-ence frame of the n -body integration. The lack ofprecision with the expansion in absolute inclina-tion was expected since absolute inclinations arenot small as ∆ I increases. However, it is interest-ing to see that even at ∆ I = 60 ◦ , the Hamiltonianexpanded in absolute inclination provides reason-able evolutions. This is due to the fact that thesystem is well described by the Laplace-Lagrangeapproximation (second order in inclination and ec-centricity), and that this approximation is invari-ant by rotation. Once again, Le Verrier’s formal-ism and the vectorial approach match perfectly,and remain in very good agreement with the full n -body integration. ∆ I = ◦ time (Myr) ∆ I = ◦ ∆ I = ◦ ∆ I = ◦ time (Myr) ∆ I = ◦ n -body ∆ I = ◦ Mercury’s eccentricity Fig. 3.— Effect of a rotation by ∆ I of thewhole system on Mercury’s eccentricity simulatedwith different secular models. The latter are: (a)expansion in eccentricity and mutual inclinationin Milankovitch’s variables, (b) expansion in ec-centricity and mutual inclination using Le Ver-rier’s equations, and (c) expansion in eccentricityand absolute inclination (Laskar & Robutel 1995).Models (a) and (b) overlap and cannot be distin-guished from one another. To test our expansions in semimajor axis (29)and (33), we integrate a system composed of aplanet with mass m = 1 M J , semimajor axis a = 6 au, and initial eccentricity e = 0 . and a brown dwarf with mass m = 40 M J , semi-major axis a = 100 au, and initial eccentric-ity e = 0 . . The mass of the central star is m = 1 M (cid:12) , and the mutual inclination is initiallyset to J = 65 ◦ . With these values, the systemundergoes Lidov-Kozai oscillations whose model-ing requires the octupole order (Naoz et al. 2011).The integrations are performed with an n -bodycode, and with the secular approximations (29)and (33). In all simulations general relavity (37)is included. The results are displayed Fig. 5. Thetwo secular models are strictly equivalent, the so-lutions are thus indistinguishable. This was notthe case for the Hamiltonian expanded in mutual14 I = ◦ time (Myr) ∆ I = ◦ ∆ I = ◦ ∆ I = ◦ time (Myr) ∆ I = ◦ n -body ∆ I = ◦ Mercury’s inclination ( ◦ ) Fig. 4.— Same as Fig. 3 but with Mercury’s in-clination. In all panels, inclinations are computedwith respect to the same reference plane after per-forming a rotation of − ∆ I on the output of eachsimulation. Models (a) and (b) are overlappingand cannot be distinguished from one another.inclination and eccentricities since the quantities p ( k, (cid:96) ) and q ( k, (cid:96) ) Eq. (22), present in the Hamil-tonian (6), were themselves expanded in T , T , V , V , W , and W to get the expression (26).The comparison with the n -body integration isalso very good. The integrations differ slightlyonce the planet flips into a retrograde orbit. Thistransition does not occur at the exact same timein the n -body and the secular simulations, as aresult, the direction of the ninth kick in the evolu-tion of I is not the same in the two integrations.Nevertheless, the transition is preceded by a pas-sage through an extreme and unrealistic eccentric-ity where the periastron distance, 0.06 R (cid:12) , is muchshorter than the size of the central star. 4. Conclusion We have provided a vectorial formalism forstudying the secular motion of non-resonant con-servative gravitational systems with concentric or-bits such as planetary systems. All expressionsare analytical and expressed in terms of the vec-tors e and j , namely the Laplace-Runge-Lenz vec- secular (b)secular (a) n -bodybrown dwarf time (Myr) time (Myr) inclination ( ◦ ) eccentricity Fig. 5.— Comparison between an n -body sim-ulation and secular integrations of a hierarchicalthree-body system. The parameters are takenfrom (Naoz et al. 2011, fig. 1). The system is com-posed of a planet with mass m = 1 M J , at a = 6 au with an initial eccentricity of e = 0 . , and abrown dwarf with mass m = 40 M J , at a = 100 au and an initial eccentricity e = 0 . . The initialmutual inclination is J = 65 ◦ . The secular modelsare: (a) expansion in eccentricity and mutual incli-nation in Milankovitch’s variables and (b) expan-sion in eccentricity and mutual inclination usingLe Verrier’s equations. Note that the eccentric-ity of the planet (upper left panel) is plotted in a tanh − scale.tor and the dimensionless orbital angular momen-tum. Planet-planet interactions have been de-veloped either in eccentricity and mutual inclina-tion or in semimajor axis ratio using an algorithmadapted from Abdullah (2001). For completeness,the vectorial expression of the spin-orbit interac-tion and general relativity have also been recalled.With numerical tests, we have shown that the in-tegrations of the vectorial equations are in perfectagreement with the more standard approach rely-ing on classical elliptic elements ( a, e, I, λ, (cid:36), Ω) .The Kepler spacecraft has revealed a populationof low-inclination, low-mass, and low-period plan-ets which are not readily studied with N-bodytechniques. The present formalism can be usedto study their dynamics, including the onset ofsecular chaos. Moreover, the vectorial approachnaturally is expanded to include a companion ona hierarchical orbit. That companion could be abinary star companion, or a giant planet at a con-15iderable distance. As long as the inner planetsystem remains with coplanar orbits and low ec-centricities, our formalism remains valid. It can,in fact, be used to see whether the system will re-main stable, or whether secular chaos will lead toorbit crossings, which then must be followed withan N-body approach. Inclusion of the precessiondynamics of central star also allows us to addressthe problem of spin-orbit alignment in multiplanetsystems, which is newly observationally relevant,and which we will pursue in a follow-up paper(Boué and Fabrycky, 2014).GB thanks Philippe Robutel, Jacques Laskarand Alain Albouy for the many discussions whichhave been helpful for this study, and in particularthose about Khaled Adbulah’s PhD thesis. 16 . Explicit vectorial equations of motion In this appendix, we provide the explicit expressions of the secular equations of motion written in avectorial form. Notations are the same as in Section 2.1.2. We consider both compact and hierarchicalsystems. In each case, the parameters of the outer body are noted with a prime, while those of the innerbody are unprimed. A.1. Compact system The Hamiltonian (26) describing the secular evolution of a compact two planet system reads H = − G mm (cid:48) a (cid:48) (cid:104) c + c ( e + e (cid:48) + j · j (cid:48) − 1) + c ( e · e (cid:48) ) + c e e (cid:48) + c ( α e + e (cid:48) ) + c ( e · e (cid:48) ) + c (1 − j · j (cid:48) ) + c ( e · j (cid:48) ) + c ( j · e (cid:48) ) + c (cid:0) α (1 − j · j (cid:48) ) e + ( e · e (cid:48) ) e (cid:48) (cid:1) + c (cid:0) (1 − j · j (cid:48) ) e (cid:48) + α ( e · e (cid:48) ) e (cid:1) + c (cid:0) (1 − j · j (cid:48) )( e · e (cid:48) ) − ( e · j (cid:48) )( j · e (cid:48) ) (cid:1) (cid:105) , (A1)where the c i ’s are parameters depending on Laplace coefficients (see Tab. 3). The conservation of the orbitalangular momentum implies that Λ d j dt = − Λ (cid:48) d j (cid:48) dt (cid:48) = T , (A2)where T = − j × ∇ j H − e × ∇ e H is a torque whose expression is T = G mm (cid:48) a (cid:48) (cid:16) A j × j (cid:48) + B e × e (cid:48) + C e × j (cid:48) + D j × e (cid:48) (cid:17) . (A3)The other equations of motion, deduced from (14), are d e dt = G mm (cid:48) a (cid:48) Λ (cid:16) A e × j (cid:48) + B j × e (cid:48) + C j × j (cid:48) + D e × e (cid:48) + E j × e (cid:17) , (A4)and d e (cid:48) dt = G mm (cid:48) a (cid:48) Λ (cid:48) (cid:16) A e (cid:48) × j + B j (cid:48) × e + C e (cid:48) × e + D j (cid:48) × j + F j (cid:48) × e (cid:48) (cid:17) , (A5)with A = c − c (1 − j · j (cid:48) ) − αc e − c e (cid:48) − c ( e · e (cid:48) ) ,B = c + 2 c ( e · e (cid:48) ) + c e (cid:48) + αc e + c (1 − j · j (cid:48) ) ,C = 2 c ( e · j (cid:48) ) − c ( j · e (cid:48) ) ,D = 2 c ( j · e (cid:48) ) − c ( e · j (cid:48) ) ,E = 2 c + 2 c e (cid:48) + 4 α c e + 2 αc (1 − j · j (cid:48) ) + 2 αc ( e · e (cid:48) ) ,F = 2 c + 2 c e + 4 c e (cid:48) + 2 c ( e · e (cid:48) ) + 2 c (1 − j · j (cid:48) ) . (A6) A.2. Hierarchical system The Hamiltonian (33) describing the secular evolution of a hierarchical two planet system reads H = − G mm (cid:48) a (cid:48) (cid:104) α j (cid:48) (cid:0) j · j (cid:48) ) − (1 − e ) j (cid:48) − e · j (cid:48) ) (cid:1) + 15 α j (cid:48) (cid:16)(cid:2) (1 − e ) j (cid:48) + 35( e · j (cid:48) ) − j · j (cid:48) ) (cid:3) ( e · e (cid:48) ) − e · j (cid:48) )( j · e (cid:48) )( j · j (cid:48) ) (cid:17)(cid:105) , (A7)17here j = (cid:107) j (cid:107) = √ − e . The conservation of the orbital angular momentum implies that Λ d j dt = − Λ (cid:48) d j (cid:48) dt = T , (A8)where T = − j × ∇ j H − e × ∇ e H is a torque whose expression is T = G mm (cid:48) a (cid:48) (cid:16) A j × j (cid:48) + B e × e (cid:48) + C e × j (cid:48) + D j × e (cid:48) (cid:17) . (A9)The other equations of motion, deduced from (14), are d e dt = G mm (cid:48) a (cid:48) Λ (cid:16) A e × j (cid:48) + B j × e (cid:48) + C j × j (cid:48) + D e × e (cid:48) + E j × e (cid:17) , (A10)and d e (cid:48) dt = G mm (cid:48) a (cid:48) Λ (cid:48) (cid:16) A e (cid:48) × j + B j (cid:48) × e + C e (cid:48) × e + D j (cid:48) × j + F j (cid:48) × e (cid:48) (cid:17) , (A11)with A = 3 α j (cid:48) ( j · j (cid:48) ) − α j (cid:48) (cid:2) ( e · j (cid:48) )( j · e (cid:48) ) + ( e · e (cid:48) )( j · j (cid:48) ) (cid:3) ,B = 15 α j (cid:48) (cid:2) (1 − e ) j (cid:48) + 35( e · j (cid:48) ) − j · j (cid:48) ) (cid:3) ,C = − α j (cid:48) ( e · j (cid:48) ) + 75 α j (cid:48) (cid:2) e · j (cid:48) )( e · e (cid:48) ) − ( j · e (cid:48) )( j · j (cid:48) ) (cid:3) ,D = − α j (cid:48) ( e · j (cid:48) )( j · j (cid:48) ) ,E = 3 α j (cid:48) − α j (cid:48) ( e · e (cid:48) ) ,F = 3 α j (cid:48) (cid:2) j · j (cid:48) ) − (1 − e ) j (cid:48) − e · j (cid:48) ) (cid:3) + 75 α j (cid:48) (cid:16) (cid:2) (1 − e ) j (cid:48) + 49( e · j (cid:48) ) − j · j (cid:48) ) (cid:3) ( e · e (cid:48) ) − e · j (cid:48) )( j · e (cid:48) )( j · j (cid:48) ) (cid:17) . (A12) A.3. Spin-orbit interaction The Hamiltonian (34) governing the spin-orbit evolution of a planet orbiting an oblate star is H = Gm mJ R a (1 − e ) / (cid:0) − s · w ) (cid:1) . (A13)Once again, the conservation of the angular momentum implies that L d s dt = − Λ d j dt = T , (A14)with T = − s × ∇ s H is the torque acting on the stellar rotation. Its expression is T = 32 Gm mJ R a (1 − e ) / ( s · w ) s × w . (A15)The evolution of the eccentricity vector e , deduced from (14), is given by Λ (cid:112) − e d e dt = − Gm mJ R a (1 − e ) / (cid:20) ( s · w ) s × e + 12 (cid:0) − s · w ) (cid:1) w × e (cid:21) . (A16)18 .4. General relativity The effect of general relativity induced by the massive central star is modeled by the Hamiltonian (37) H = − µ βa c √ − e . (A17)The orbital angular momentum of the planet is conserved, and its eccentricity vector evolves according to d e dt = 3 µ βa c Λ(1 − e ) w × e . (A18)19 EFERENCES TRIP Hansen, P. A. 1853, Abhandl. d. K. S. Ges. d.Wissensch., IV, 182Hill, G. W. 1875, The Analyst, 2, 161Holman, M., Touma, J., & Tremaine, S. 1997, Na-ture, 386, 254 Innanen, K. A., Zheng, J. Q., Mikkola, S., & Val-tonen, M. J. 1997, AJ, 113, 1915Kaib, N. A., Raymond, S. N., & Duncan, M. J.2011, ApJ, 742, L24Katz, B., Dong, S., & Malhotra, R. 2011, PhysicalReview Letters, 107, 181101Kozai, Y. 1962, AJ, 67, 591Kratter, K. M., & Perets, H. B. 2012, ApJ, 753,91Lai, D. 2014, MNRAS, 440, 3532Lambeck, K. 1988, Geophysical geodesy : the slowdeformations of the earth (Oxford [England] :Clarendon Press ; New York : Oxford Univer-sity Press, 1988.)Laskar, J. 1984, Thèse, Observatoire de Paris http://tel.archives-ouvertes.fr/tel-00702723 —. 1989, Nature, 338, 237—. 1990, Icarus, 88, 266—. 2008, Icarus, 196, 1Laskar, J., & Boué, G. 2010, A&A, 522, A60Laskar, J., & Robutel, P. 1995, Celestial Mechan-ics and Dynamical Astronomy, 62, 193Le Verrier, U. J. J. 1855, Annales de l’Observatoireimpérial de Paris: mémoires (Impr. Gauthier-Villars)Li, G., Naoz, S., Kocsis, B., & Loeb, A. 2014, ApJ,785, 116Libert, A.-S., & Delsate, N. 2012, MNRAS, 422,2725Libert, A.-S., & Henrard, J. 2007, Icarus, 191, 469—. 2008, Celestial Mechanics and Dynamical As-tronomy, 100, 209Libert, A.-S., & Tsiganis, K. 2009, MNRAS, 400,1373—. 2011, MNRAS, 412, 2353Lidov, M. L. 1962, Planet. Space Sci., 9, 71920ithwick, Y., & Naoz, S. 2011, ApJ, 742, 94Lithwick, Y., & Wu, Y. 2011, ApJ, 739, 31Mardling, R. A. 2010, MNRAS, 407, 1048Michtchenko, T. A., Ferraz-Mello, S., & Beaugé,C. 2006, Icarus, 181, 555Migaszewski, C., & Goździewski, K. 2009, MN-RAS, 395, 1777Milankovitch, M. 1939, Bull. Serb. Acad. Math.Nat. A, 6, 1Morbidelli, A., Tsiganis, K., Batygin, K., Crida,A., & Gomes, R. 2012, Icarus, 219, 737Mugrauer, M., Neuhäuser, R., Mazeh, T., et al.2006, Astronomische Nachrichten, 327, 321Nagasawa, M., & Ida, S. 2011, ApJ, 742, 72Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A.,& Teyssandier, J. 2011, Nature, 473, 187—. 2013a, MNRAS, 431, 2155Naoz, S., Farr, W. M., & Rasio, F. A. 2012, ApJ,754, L36Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013b,ApJ, 773, 187Petrovich, C. 2014, ArXiv e-prints,arXiv:1405.0280Raghavan, D., McAlister, H. A., Henry, T. J.,et al. 2010, ApJS, 190, 1Rosengren, A. J., & Scheeres, D. J. 2014, CelestialMechanics and Dynamical Astronomy, 118, 197Saleh, L. A., & Rasio, F. A. 2009, ApJ, 694, 1566Souriau, J. M. 1969, Structure des systèmes dy-namiques (Dunod, Paris)Takeda, G., Kita, R., & Rasio, F. A. 2008, ApJ,683, 1063Terquem, C., & Papaloizou, J. C. B. 2002, MN-RAS, 332, L39Tisserand, F. 1889, Traité de Mécanique Céleste(Paris: Gauthier-Villars) Touma, J. R., Tremaine, S., & Kazandjian, M. V.2009, MNRAS, 394, 1085Tremaine, S., & Dong, S. 2012, AJ, 143, 94Tremaine, S., Touma, J., & Namouni, F. 2009, AJ,137, 3706Tremaine, S., & Yavetz, T. 2013, ArXiv e-prints,arXiv:1309.5244Triaud, A. H. M. J. 2011, A&A, 534, L6Triaud, A. H. M. J., Collier Cameron, A., Queloz,D., et al. 2010, A&A, 524, A25Udry, S., & Santos, N. C. 2007, ARA&A, 45, 397Veras, D., & Tout, C. A. 2012, MNRAS, 422, 1648Winn, J. N., Fabrycky, D., Albrecht, S., & John-son, J. A. 2010, ApJ, 718, L145Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109—. 2013, ApJ, 772, 74Wu, Y., & Murray, N. 2003, ApJ, 589, 605Wu, Y., Murray, N. W., & Ramsahai, J. M. 2007,ApJ, 670, 820Yoder, C. F. 1995, in Global Earth Physics: AHandbook of Physical Constants, ed. T. J.Ahrens, 1Zucker, S., & Mazeh, T. 2002, ApJ, 568, L113