Secular dynamics in extrasolar systems with two planets in mutually inclined orbits

Abstract We revisit the problem of the secular dynamics in two-planet systems in which the planetary orbits exhibit a high value of the mutual inclination. We propose a ‘basic hamiltonian model’ for secular dynamics, parameterized in terms of the system’s Angular Momentum Deficit (AMD). The secular Hamiltonian can be obtained in closed form, using multipole expansions in powers of the distance ratio between the planets, or in the usual Laplace-Lagrange form. The main features of the phase space (number and stability of periodic orbits, bifurcations from the main apsidal corotation resonances, Kozai resonance etc.) can all be recovered by choosing the corresponding terms in the ‘basic Hamiltonian’. Applications include the semi-analytical determination of the actual orbital state of the system using Hamiltonian normalization techniques. An example is discussed referring to the system of two outermost planets of the ν-Andromedae system.


Introduction
The dynamics of exoplanetary systems with two planets in orbits with non-zero mutual inclination is a very interesting topic in view of the discovery, in the last 20 years, of several such systems (see Naoz (2016) and references therein, as well as Libert & Tsiganis (2009)). In the present short review we do a preliminary discussion of the phase space structures observed in such systems; a more complete study will be presented in a separate work. In particular we discuss below two distinct characteristic regimes. The first, called the nearly-planar regime, is characterized by small values of the mutual inclination, high values of the planetary eccentricities and by the predominance of orbits (periodic or quasiperiodic) linked to apsidal corotations. The second regime, called Kozai-Lidov regime, is characterized, instead, by a high value of the mutual inclination, small values of the eccentricities and by the instability of the circular orbit for the inner planet. We will finally discuss the sequence of bifurcations that connect the two regimes.

Hamiltonian model
The Hamiltonian of the three-body problem in Poincaré heliocentric canonical variables takes the form: Keplerian part (2.1) where m 0 = mass of the star, and m i , p i , r i , i = 2, 3 are the masses, barycentric momenta and heliocentric position vectors of the planets. Starting from (2.1), a secular Hamiltonian is arrived at by averaging the above Hamiltonian with respect to all short period terms. We consider two distinct types of expansion, i.e. i) series expansions in powers of the planets' eccentricities and inclinations and ii) closed-form averaging. In both cases, we end up with a secular Hamiltonian model of the form In the Hamiltonian (2.2) the 'fast' angles λ i = M i + i , i = 2, 3 , are ignorable, a fact implying the constancy of the semi-major axes under the secular model. The total angular momentum L = r 2 × p 2 + r 3 × p 3 is an exact first integral of the system, a fact implying that the Hamiltonians (2.1) and (2.2) depend on the angles Ω 2 , Ω 3 only through the difference Ω 2 − Ω 3 . Using modified Delaunay canonical variables the secular Hamiltonian H sec (Γ 2 , Γ 3 , Θ 2 , Θ 3 , γ 2 , γ 3 , ϑ 2 , ϑ 3 ) has 4 degrees of freedom. However, the existence of two independent integrals in involution (i.e. the components L z and L plane of the total angular momentum L ) allows to reduce the number of degrees of freedom by two, a process known as 'Jacobi's reduction of the nodes'. We propose a novel method to perform Jacobi's reduction with respect to the one followed by Libert & Henrard (2007). Our method leads to an explicit analytic control of all small parameters appearing in the problem, by taking advantage of the symmetries of the Hamiltonian expressed in Keplerian heliocetric elements with respect to the Laplace reference frame. In particular setting, as usual, Ω 3 − Ω 2 = π and observing that the Hamiltonian depends only on the mutual inclination, H sec = H sec (a 2 , a 3 , e 2 , e 3 , cos(i 2 + i 3 ), ω 2 , ω 3 , Ω 3 − Ω 2 = π ), we introduce two book-keeping parameters, ε and η (with numerical values ε = η = 1 ), via the relations In view of these book-keeping parameters, the Hamiltonian takes the form where, in the last passage, all powers η s1 with s 1 1 are replaced simply by η (this is a possible simbolic operation since η = 1 ). Finally, setting where C = L z and G j = m j G m 0 a j 1 − e 2 j j = 2, 3 , the Hamiltonian takes the form H sec (a 2 , a 3 , e 2 , e 3 , ω 2 , ω 3 ; We note that Eq. (2.4) explicitly keeps track of all the small quantities of the problem: these are the planetary eccentrities e 2 , e 3 and the inclinations i 2 , i 3 . However, in our approach the dependence of the Hamiltonian on the inclinations is accounted for through powers of the quantities sin(i 2 ) sin(i 3 ) and cos(i 2 ) cos(i 3 ) − 1 . Finally it is possible to express the secular Hamiltonian in the form of a polynomial H sec (X 2 , X 3 , Y 2 , Y 3 ; AMD), where the canonical pairs (X 2 , Y 2 ), (X 3 , Y 3 ) are Poincaré variables We checked the precision of each type of expansion against numerical experiments; in particular we compared the Poincaré section Y 3 = 0 ,Ẏ 3 0 obtained from the Hamiltonian expanded as in both methods above. This allows to specify which is a sufficient order of truncation of the multipolar expansion and of the small parameters; Figure 1 shows an example of such comparison, allowing to see that we need to reach degree 5 of the multipolar expansion to obtain the same phase portraits as in the Laplace expansion up to order 10 in the eccentricities.

Phase portraits
In order to analyze the most important phenomena in the phase portrait of the Hamiltonian H sec (X 2 , X 3 , Y 2 , Y 3 ; AMD) , we consider, as a representative example, the value AMD = 0.016 in a fictitious system with two planets same as the planets c and d of the system ν-Andromedae. In particular we chose m 0 = 1.31 M , m 2 = 13.98 M J , m 3 = 10.25 M J , a 2 = 0.829 AU , a 3 = 2.53 AU (according to Table 13 of McArthur et al. (2010)) and e 2 = 0.2445 , e 3 = 0.316 , i 2 = 11.347 • and i 3 = 25.609 • (according to Table 1 of Deitrick et al. (2015)). Figure 2 shows the Poincaré section Y 3 = 0 ,Ẏ 3 ≥ 0 for different values of the energy E = H sec . Since the AMD is fixed, altering the values of the planets' eccentricities implies that the inclinations also change to keep the AMD constant to its pre-selected value. The maximum values of e 2 and e 3 allowed for a specific value of the AMD can be computed using the Lagrange multiplier method. In the same way we find the allowed region (X 2min , X 2max ) for the phase portrait corresponding to each section.

Nearly-planar regime
The first two panels in the Figure 2 are representative of the phase portrait obtained in the nearly-planar regime. This is similar to the phase portrait found in the exact planar case, in which the averaged Hamiltonian turns to be integrable (see Figure 3c). In fact the 3D Hamiltonian can be decomposed as Of particular importance in the integrable part H int are the periodic orbits called antialigned (mode A ) and aligned (mode B ) apsidal corotation (see Laughlin et al. (2002), Lee & Peale (2003), Beaugé et al. (2003)). In order to compute all possible (symmetric or asymmetric) apsidal corotations for given energy and angular momentum it is particularly convenient to express the integrable part of the Hamiltonian H int in Weyl variables, defined as where X i , Y i i = 2, 3 are the canonical Poincaré variables defined above. These variables satisfy the Poisson algebra relations {σ i , σ j } = −2ε ijk σ k , where ε ijk is the Levi-Civita Note the absence of a separatrix between the two stable modes. This is due to the fact that the system's reduced phase-space has the topology of the 3D-sphere.
symbol, when i, j, k = 1, 2, 3 and ε ijk = 0 if one of the i, j, k is equal to 0 . We easily verify that the Hamiltonian H int does not depend on σ 2 . Consider the surfaces (3.1) It can be proved that the periodic orbits of H int are given by the points of tangency of the surfaces S σ0 and C σ0,E . Using this property we can compute analytically all the apsidal corotations corresponding to a fixed energy E (varying the angular momenta of the planets, i.e. σ 0 ), or fixed σ 0 (varying E ). Figure 3 shows a particular example.
Having specified, now, the coordinates of the apsidal corotations, it is possible to apply a Birkhoff normal form of the Hamiltonian in order to compute the quasi-periodic orbits around one of the two modes. Details of this construction are given in a separate study. We emphasize the method's utility in order to provide a semi-analitycal representation of the long term time series of the orbital elements for the planetary trajectories.

Sequences of bifurcations
As shown in Figure 2, the nearly-planar regime is followed by a sequence of bifurcations occurring for a fixed value of the AMD as the energy increases (from left to right), implying that the maximum allowed mutual inclination between the planets also increases. Starting with the basic apsidal corotations (A, B), a saddle-mode bifurcation generates the orbits C 1 , C 2 which correspond to an orbital configuration with nonzero mutual inclination. Furthermore, as the mutual inclination increases, the orbit C 2 becomes unstable by the "Kozai mechanism"(see Kozai (1962)), as shown in the last picture of Figure 2.

Kozai-Lidov regime
In the last two pictures of Figure 2 we can see the transition of C 2 from stable to unstable, accompanied by a large volume of the trajectories around C 2 becoming chaotic. The transition occurs at critical values of L z , or, equivalently, of the mutual inclinations; for example the limiting inclination in the case of an inner test particle and an outer planet in circular orbits is cos −1 3 5 ∼ 39 • .2 (see Naoz (2016) for a review). These limits