The phase-space architecture in extrasolar systems with two planets in orbits of high mutual inclination
R<EMAIL_ADDRESS>and C<EMAIL_ADDRESS>
Dipartimento di Matematica Tullio Levi-Civita, Università degli
Studi di Padova

via Trieste 63, 35121 Padova

We revisit the secular 3D planetary three-body problem with the aim to provide a unified formalism allowing for a convenient parametric representation of all basic phenomena sculpting the form of the phase space as the mutual inclination between the two planets' orbits is set to gradually larger values, starting from zero (planar case) and ending to a value beyond the limit of the Lidov-Kozai instability. Our main findings are: i) on methodological ground, we propose an algebraic `book-keeping' technique, incorporated to the treatment of the secular Hamiltonian already at the level of the Jacobi reduction of the nodes, which allows to decompose the secular Hamiltonian in the form $\Hscr_{sec}=\Hscr_{planar}+\Hscr_{space}$, with $\Hscr_{space}$ collecting all terms depending on the planets mutual inclination $\i_{mut}$. This technique can be implemented both in expansions of the Hamiltonian in multipoles (Legendre polynomials) or in power of the planetary eccentricities (Laplace-Lagrange). We numerically compare several such models, aiming to establish an appropriate value of the minimum multipole truncation order $N_P$, as well as of the minimum book-keeping truncation order $N_{bk}$ (accounting simultaneously for powers of the planetary eccentricities, inclinations, as well as the system's Angular Momentum Deficit ($\AMD$)), such that the form of obtained phase portrait stabilizes for truncations at orders beyond the minimum ones. This analysis is made in a regime far from any hierarchical limit either in the planetary mass ratio or distance ratio. ii) On dynamical ground, we explore the transition, as the level of mutual inclination increases, from a `planar-like' to a `Lidov-Kozai' regime, both characterized by the kind of periodic orbits which dominate the structure of the phase space at the corresponding values of the energy. Using a typical case of non-hierarchical parameters (pertinent to the well known $\upsilon$-Andromedae system), we find that the structure of the phase portraits of the integrable secular dynamics of the planar case is reproduced to a large extent and for a wide range of energies also in the 3D case. We provide a semi-analytical criterion allowing to estimate the extent in energies and level of mutual inclination up to which dynamics remains nearly-integrable. In this regime, we propose a normal form method by which the basic periodic orbits of the nearly-integrable regime (modes A and B), can be computed semi-analytically. These modes generalize the apsidal corotation orbits (A: anti-aligned, B: aligned) of the planar case. On the other hand, we show how, as the energy increases, the system gradually moves from the nearly-integrable `planar-like' to the `Lidov-Kozai' regime. The latter is dominated by two different families of nearly circular and highly inclined periodic orbits ($C_1$ and $C_2$), of which $C_2$ becomes unstable via the usual Lidov-Kozai mechanism. We demonstrate how the $C_1$ and $C_2$ families are connected to the $A$ and $B$ families via a chain of saddle-node and pitchfork bifurcations. We provide an analytical criterion for the destabilization of the family $C_2$ obtained in terms of the system's orbital parameters (value of the eccentricity of the outer planet as well as value of the $\AMD$), obtained in the framework of the quadrupolar approximation. Finally, we perform a numerical study on the form of the phase portraits for different mass ratios as well as semi-major axis ratios of the two planets, aiming to establish how generic are the phenomena reported above as the systems parameters are chosen close to one or more hierarchical limits.

§ INTRODUCTION

In the present paper we revisit the planetary three body problem in Poincaré heliocentric canonical variables governed by the Hamiltonian

\begin{equation}
\label{ham.3BP}
\Hscr(\vet{r_2},\vet{r_3},\vet{p_2},\vet{p_3})=\frac{\vet{p_2}^2}{2\, m_2}-\frac{\mathcal{G}\, m_0\, m_2}{r_2}+\frac{\vet{p_3}^2}{2 \,m_3}-\frac{\mathcal{G}\, m_0\, m_3}{r_3}+\frac{(\vet{p_2}+\vet{p_3})^2}{2\, m_0}-\frac{\mathcal{G}\, m_2\, m_3}{|\vet{r_2}-\vet{r_3}|}\, ,
\end{equation}

where $m_0 \gg m_2,\, m_3\, $. Thus $m_0$ represents the mass of a star, and $m_{i}\,$, $\vet{p_i}\,$, $\vet{r_{i}}\,$, $i=2,3$ the masses, barycentric momenta and heliocentric position vectors of two planets orbiting the star, at distances $r_i=|\vet{r_i}|\,$.

Our focus below will be on a systematic study of how the structure of the phase space of the above Hamiltonian system is altered as the planetary masses and distances and, most importantly, the mutual inclination between the two planets' trajectories is varied. Answering this question in the framework of the three-body problem is a key step towards understanding the phase-space architecture in planetary systems with two or more planets in orbits with high mutual inclination. A well-studied case of the latter is the $\upsilon$-Andromedae system. This is a double star system with four planets (b, c, d and e) orbiting one of the stars. Since the mass $m_b$ is much smaller than $m_c\, , m_d$ the motion of the innermost planet b can be modelel to a good approximation via a restricted four-body problem, with planets c, d providing the main perturbations (this motivates the use of the indices 2,3 in the Hamiltonian (<ref>) for the masses and canonical variables of the planets); the long-term motion of the innermost planet $(m_1,\vet{r}_1,\vet{p}_1)$ is considered in a separate work (Mastroianni and Locatelli, in preparation). More detailed models can include planet e as well as the second star. At any rate, all the above models require providing first a good analytical model for the orbits of the giant planets of the system, say, for example, the planets c and d, whose masses, as estimated by observations, are larger than $10\,M_{J}\,$ (see e.g [Deitrick et al., 2015] and [McArthur et al., 2010]).

Now, the available data for the orbital parameters of an extrasolar system are typically affected by wide observational error bars. As will be discussed in detail below, one important problem with the uncertainties in the observations stems from the fact that even small changes in a system's estimated parameters, consistent with the observations, may imply a drastic change in the type of orbital state in which the observed system is assumed to have been settled. We will argue below that this sensitive dependence on  available parameter estimates affects mostly those predictions referring to the secular (long-term) evolution of the orbital state, i.e., the variations of the planets' eccentricity and inclination vectors which take place in timescales of the order of $\sim 10^2-10^4$ orbital periods. Notwithstanding these long timescales, characterizing the whole variety of possible stable secular orbital states of exoplanetary systems can be  relevant also in the interpretation of short-in-time observations: most importantly, it can serve the purpose of constraining observational uncertainties on the basis of stability considerations, i.e., indicating which subdomains in parameter space favor the long-term stability of the planetary orbits [Petit et al., 2018].

Our study in the present paper focuses on one phenomenon, whose role appears central to the aim of classifying and characterizing the variety of possible secular orbital states in 3D planetary systems: this is the chain of bifurcations of periodic orbits which mark the transition, as the planets' mutual inclination increases, from the apsidal corotation to the Lidov-Kozai regime.

As is well known, the apsidal corotation (AC) states (see, for example, [Laughlin et al., 2002], [Beaugé et al., 2003], [Lee and Peale, 2003]) provide the backbone of the phase space of secular motions for two planets in coplanar orbits. Formally, the apsidal corotations are the continuation, in the nonlinear regime, of the linear (Laplace-Lagrange) normal mode solutions of the planar two-planet secular Hamiltonian. Physically, they represent periodic orbits along which the pericenters of the two planets constantly precess by remaining either always anti-aligned (state A) or always aligned (state B). Periodic orbits having the same property in the exact (non-averaged with respect to fast angles) Hamiltonian can also be found [Hadjidemetriou, 1975], [Hadjidemetriou, 2006] once the dynamics is regarded in a frame rotating with angular velocity equal to the (common) precession rate of the apsides.

Given their importance in the planar case, a natural question regards how the family of apsidal corotation orbits is continued, as well as what is the dynamical role played by the AC periodic, and nearby quasi-periodic, orbits, when we pass from the coplanar to the 3D planetary orbital configuration. Both the above questions have been addressed in the literature using diverse formalisms (see [Beaugé et al., 2012] and references therein, as well as the discussion below). Owing to reasons explained in detail in section <ref> below, in our present study we employ our own-proposed formalism, in which, by applying a so-called book-keeping technique already at the level of Jacobi's reduction of the nodes [Jacobi, 1842]), we arrive at a natural decomposition of the 3D secular Hamiltonian of a system with fixed Angular Momentum Deficit (AMD; see [Laskar and Robutel, 1995]) as

\begin{equation}\label{h3ddecompo}
\Hscr_{sec}=\Hscr_{planar}(\mathbf{X},\mathbf{Y})+\Hscr_{space}(\mathbf{X},\mathbf{Y};\AMD)~~.
\end{equation}

In Eq.(<ref>), $(\mathbf{X},\mathbf{Y})$ are Poincaré canonical variables for the two planets ($(X_2,Y_2)$, $(X_3,Y_3)$ are approximately proportional to the planets' eccentricity vectors). The Angular Momentum Deficit is defined by

\AMD=L_2+L_3-L_z

where $L_2$, $L_3$ are the angular momenta of the circular orbits at semi-major axes $a_2,a_3$ equal to those of the two planets, and $L_z$ the modulus of the total angular momentum normal to the system's Laplace plane.

The methodological benefits from working with a decomposition of the secular Hamiltonian as in Eq.(<ref>) stem from the following properties (see sections 2 and 3 for details):

i) $\Hscr_{planar}(\mathbf{X},\mathbf{Y})$ is an integrable Hamiltonian, the quantity $J=(X_2^2+Y_2^2+X_3^2+Y_3^2)/4$ being a second integral independent of the energy.

ii) For every permissible value of the energy $\Escr$, all the orbits $(\mathbf{X}(t),\mathbf{Y}(t))$ under Hamilton's equations with the Hamiltonian $\Hscr_{planar}$ are confined to the Laplace plane, i.e., they have zero mutual inclination at the given (and fixed in advance) level of AMD.

iii) The Hamiltonian $\Hscr_{planar}$ can be explicitly shown to admit two periodic orbits (called below the `modes' A and B) which correspond to the anti-aligned and aligned apsidal corotation states for the two planets in the planar case.

iv) The Hamiltonian $\Hscr_{space}$ can be further decomposed as $\Hscr_{space}=\Hscr_{0,space}+\Hscr_{1,space}$, where $\Hscr_{0,space}$ also admits $J=(X_2^2+Y_2^2+X_3^2+Y_3^2)/4$ as a second integral. Then, the Hamiltonian

\begin{equation}\label{def.h0}
\Hscr_{int}=\Hscr_{planar}+\Hscr_{0,space}
\end{equation}

is integrable, and it has a similar formal structure as the Hamiltonian $\Hscr_{planar}$. In particular, the existence of periodic orbits of the type A and B can be explicitly demonstrated for the Hamiltonian $\Hscr_{int}$ using the integrability property, in the same way as for the Hamiltonian $\Hscr_{planar}$.

v) Focusing, now, on the full secular Hamiltonian

\begin{equation}\label{hamful}
\Hscr_{sec}=\Hscr_{int}+\Hscr_{1,space}
\end{equation}

a Birkhoff-like normal form construction (section <ref>) allows to demonstrate (up to an exponentially small error) that the periodic orbits A,B (accompanied by neighboring quasi-periodic orbits) continue to exist in the full 3D regime. Hence, we call these orbits the 3D apsidal corotation states.

vi) We parameterize the level of non-coplanarity (mutual inclination $\i_2+\i_3$) of the system in terms of the (constant) energy level of the full Hamiltonian $\Escr=\Hscr_{sec}$. One way to regard the connection between the value of the energy $\Hscr_{sec}$ and the level of non-coplanarity of the orbits is by noting that the energy grows in absolute value nearly as a quadratic function of the planetary eccentricities, i.e., nearly proportionally to a linear combination of $\e_2^2$ and $\e_3^2$. Thus, restricted to the chosen Poincaré surface of section (see below), the constant energy condition $\Escr=\Hscr_{sec}$ yields an ellipsoid-like surface. Also (see section <ref>), for fixed $\AMD$, the condition of constant mutual inclination yields the contour of a function also quadratic in the orbital eccentricities. Then, as shown in detail in section <ref>, for any fixed value of the energy $\Escr$, there are two values of the mutual inclination $\i_{mut}=\i_2+\i_3$, namely $\i_{mut}^{min}(\Escr)$, $i_{mut}^{max}(\Escr)$ such that, the corresponding contours of constant mutual inclination come tangent to the ellipsoidal surface of constant energy $\Escr=\Hscr_{sec}$. As a consequence, for any value of the mutual inclination in the interval $\i_{mut}^{min}(\Escr)\leq \i_{mut}\leq \i_{mut}^{max}(\Escr)$ there is an energetically permissible domain of motions. By the geometric properties of the above considerations, we find that both $\i_{mut}^{min}(\Escr)$ and $\i_{mut}^{max}(\Escr)$ are monotonically increasing functions of the energy. Thus, selecting a particular level value of the energy $\Escr$ fixes the overall level of mutual inclinations allowed at the energy $\Escr$. The minimum possible energy $\Escr_{min}$ is defined by the condition that the level surface $\Hscr_{sec}=\Escr_{min}$ is tangent to the level surface of minimum possible mutual inclination $\i_{mut}=0$ (planar case), while all other points of the surface $\i_{mut}=0\,$, except for the points of tangency, are in the interior of the surface $\Hscr=\Escr_{min}$. Solving for these conditions allows to specify the value of $\Escr_{min}$.

In conclusion, the level of mutual inclination for all orbits increases in general with the quantity $\delta \Escr=\Escr-\Escr_{min}$, where $\Escr_{min}\leq \Escr\leq 0$. Thus, the lowermost value of $\delta\Escr$ is $\delta\Escr_{min}=0$, while the highermost limit is $\delta\Escr_{max}=-\Escr_{min}$, i.e., $\Escr=0$. At this latter limit, the available phase space shrinks to a point with $\e_2(t)=\e_3(t)=0$. This means a unique possible orbital configuration of two mutually-inclined circular planetary orbits. We call this the limiting trajectory of the Lidov-Kozai regime. Actually, we will see in section <ref> that, for negative energies $\Escr$ close to zero, the phase space acquires a structure reminiscent to the one of the non-integrable Lidov-Kozai case of the restricted three-body problem (i.e. non intersecting trajectories examined in a higher than quadrupolar development of the disturbing function, see [Libert et al., 2011]). The typical behavior of the trajectories in the Lidov-Kozai regime is to (quasi-)periodically exchange eccentricity for mutual inclination. However, we will see that such a coarse illustration of the dynamics is rather simplistic in the case of non-hierarchical (in distance ratios or masses) two-planet systems; in reality, the dynamics around the central Lidov-Kozai periodic orbit is highly unstable and the corresponding phase space turns to exhibit strong chaos.

In summary, our focus in the present paper will be on describing, with sufficient detail, the observed transition in the structure of the phase space, as the parameter $\delta \Escr$ increases from its lowermost limit, corresponding to a nearly planar orbital configuration, to the highermost limit, corresponding to nearly circular orbits with a high degree of mutual inclination.

As expected in the study of any dynamical system, structural changes in the phase space are associated with the birth, bifurcations and stability evolution of the most important periodic orbits of the system. We already mentioned that these are the ACs, in the nearly planar limit, and the Lidov-Kozai orbits, in the maximum mutual inclination limit. We follow the evolution and the connections between these periodic orbits, as well as other ones emerging in half the way between the two limiting regimes. In our study we work with analytical estimates, as well as with a numerical example inspired by the $\upsilon$-Andromedae system. This is an example of system being reasonably far from any hierarchical limit: the estimated mass ratio of the two planets is $m_2/m_3\simeq 1.3$, while the estimated semi-major axis ratio is $a_2/a_3\simeq 0.3$. On the other hand, section 4 contains results from a numerical investigation referring to different choices for the mass and semi-major axis ratios, representing every one of all possible cases of hierarchical models that could arise in the problem under consideration.
As a rough guide, the phenomena discussed below should cover most cases of interest in the range of mass ratios $1/10\leq m_2/m_3 \leq 10\,$, distance ratios $ 1/7\leq r_2/r_3\leq 1/3 $ and mutual inclination $0\leq \i_{mut}\leq 45^\circ\, $.

Few more words about our hereby proposed formalism and method of study: the spectrum of methods and formalisms proposed for the study of the secular planetary three body problem is nearly as wide as the literature itself on the subject. Different proposals distinguish between cases in which the three body problem dealt with is considered hierarchical (e.g. in the distances $r_2/r_3 \ll 1\,$; see [Brouwer, 1959], [Ford et al., 2000], [Migaszewski and Goździewski, 2011]), or non-hierarchical (see [Henrard and Libert, 2004], [Naoz et al., 2013]). The choice of variables and/or proposed representation of the phase space of the system has often been motivated by whether the focus of a particular study is on phenomena related to apsidal-corotations (see [Beaugé et al., 2003], [Laughlin et al., 2002], [Lee and Peale, 2003]), or the Kozai instability (see [Kozai, 1962], [Libert and Delsate, 2012], [Libert and Tsiganis, 2009], [Lithwick and Naoz, 2011]).
In addition, different formalisms stem from choices related to:

i) the type of Hamiltonian expansion: this can be performed as the usual (Laplace-Lagrange) series expansion in powers of the planets' eccentricities and inclinations (see for example [Murray and Dermott, 1999] and [Libert and Henrard, 2007]) or as a Legendre multipolar expansion (see [Naoz, 2016], [Libert and Delsate, 2012] and references therein). One can note that, although the two types of expansions are equivalent in the limit of infinite order of the expansion, different truncations (even with the same type of series) can lead to quantitatively, and even qualitatively, different results. Libert & Henrard in [Libert and Henrard, 2007] have discussed the question of the correct order of truncation in the framework of the Laplace-Lagrange expansion. Mogaszewki & Goździewki and Naoz (see [Migaszewski and Goździewski, 2011], [Naoz et al., 2013] and [Naoz, 2016]) discuss various truncated multipolar models for spatial hierarchical systems. To our knowledge, there is no literature on a comparison between the results obtained by truncated models with the two types of expansion for non-hierarchial spatial systems.

ii) Method of averaging: different secular models are obtained by a different choice of method for averaging the Hamiltonian with respect to the fast angles. Such methods include: a) averaging “by scissors”,  i.e., by just dropping-off the Hamiltonian all fast periodic terms (see, for example, [Libert and Henrard, 2007], [Henrard and Libert, 2004]). b) “Closed-form” averaging (e.g. [Migaszewski and Goździewski, 2011]). This method has the benefit of avoiding expansions in the orbital eccentricities (whose convergence becomes limited due to the limit in the series inversion of Kepler's equation, see, for example [Szeto and Lambeck, 1982]). However, closed-form averaging can only be performed after a multipolar expansion of the Hamiltonian. Thus the method is particularly suited for systems hierachical with respect to the planetary distances, while its precision in the case of non-hierarchical systems is an open issue (see some results in section <ref>). c) Numerical computation (e.g. by Gauss' method) of the quadratures involved in the averaging (see [Thomas and Morbidelli, 1996], [Migaszewski and Goździewski, 2008]). This method is particularly suited for systems with intersecting trajectories, owing to theorems (see, for example, [Gronchi and Milani, 1998]) establishing the continuity of the secular equations of motion at the points of intersection, which are singularities of the integrand functions appearing in the quadratures. d) Elimination of the mean anomalies from the Hamiltonian via a canonical transformation (see [Brouwer, 1959] or [Naoz et al., 2013]). From the theoretical point of view, the use of a canonical transformation is imperative when precision of second order in the planetary masses is sought for in the secular model. It should be stressed that while $\Oscr(m^2)$ expansions are straightforward to obtain in the framework of the Laplace-Lagrange series (see [Locatelli and Giorgilli, 2000]), their counterpart in the form of closed-form series is an open question. In fact, closed-form averaging in the framework of the three-body problem requires the use of some `relegation' technique (see [Palacián et al., 2006]), or alternatives as those recently proposed in [Cavallari and Efthymiopoulos, 2022]).

iii) Choice of coordinates. Several sets of variables have been proposed as convenient for the  visualization and study of the phase space of secular motions. Examples are:
$(\e_2\,\sin(\omega_2),\,\e_3\,\sin(\omega_3))$ (see for istance [Libert and Henrard, 2007], [Migaszewski and Goździewski, 2011], [Libert and Tsiganis, 2009]);
$(\e_j\cos(\Delta\varpi),\,\e_j\sin(\Delta\varpi))\,$, $j=2,\,3\,$ (see [Michtchenko and Malhotra, 2004]);
$(\e_2\,\cos(2\omega_2),\,\e_3\cos(\Delta\omega))$ (see [Libert and Delsate, 2012]);
$(\e_2\,\cos(\Delta\varpi),\,\e_3\,\cos(2\,\omega_2)\,)$ (see [Migaszewski and Goździewski, 2011], [Michtchenko et al., 2006]);
$(\e_j\,\cos(\omega_j),\,\e_j\,\sin(\omega_j))\,$ (see [Thomas and Morbidelli, 1996]).
Among the motivations behind the choice of a particular set of variables are: a) the treatment of singular cases (e.g. $\e_2=0$, $\e_3=0$, in which the longitude of the perihelium is no longer defined), b) the possibility to include all the main families of possible orbital states in one plot. As a characteristic example of the latter, see figure 3 of [Michtchenko et al., 2006]; one can remark, however, the complexity involved in properly deciphering the information given in that figure, which is evident from the accompanying caption.

Our own choice regarding the points (i) to (iii) above is: (i) and (ii) we base most of our results on a closed-form averaging of a multipolar expansion (in powers of the distance ratio $r_2/r_3<1$) of the Hamiltonian (<ref>). Besides its compact form, convenient for analytic studies (as, for example, in subsection <ref>), such a model, truncated at a sufficiently high multipole order, circumvents the problem of slow convergence of the Laplace-Lagrange series for highly eccentric orbits, without compromising precision even far from the hierarchical limit (e.g. for $r_2/r_3\sim 0.3\,$). In order to specify the suitable order of multipole  truncation, in section <ref> we make a comparison of the phase portraits obtained via the secular Hamiltonian  arrived at by the above method versus those obtained by a scissors' averaging of the Laplace-Lagrange series truncated at order 10 in the eccentricities. Also, in obtaining the final Hamiltonian we introduce a book-keeping procedure for the Jacobi reduction of the nodes, which, as mentioned already, leads to a convenient decomposition of the Hamiltonian as in Eq.(<ref>). (iii) We illustrate all phase portraits using the usual Poincaré surface of section $(\e_2\cos\omega_2,\e_2\sin\omega_2)$ every time when $\omega_3=\pi$, $\dot{\omega_3}\geq 0$. The sequence of canonical transformations leading to such a representation of phase portraits is described in section <ref>. Owing to the conservation of angular momentum, one can easily see ([Cushman and Bates, 1997]) that the phase space of the integrable models $\Hscr_{planar}$ and $\Hscr_{int}$ is the sphere $S^2$ (instead of $\mathbb{R}^4$, as generically true for Hamiltonian systems with two degrees of freedom). Furthermore, as analyzed in section <ref>, some points of the sphere `inflate' to curves in the Poincaré surface of section defined as above. Finally, using an appropriate set of good variables whose Poisson algebra admits the angular momentum as a Casimir, the two modes A and B are seen to be separated by a meridian circle in the sphere $S^2$, which however does not correspond to a dynamical separatrix (since the integrable model contains no unstable periodic orbits). We devote some effort to carefully describe these phenomena, which are often found to generate some confusion in the literature when use is made of some of the components of the planets' eccentricity vectors to describe the structure of the problem's phase space.

The paper is organized as follows. Section <ref> describes all the steps followed in order to arrive at the finally adopted secular Hamiltonian model. This includes the introduction of an appropriate `book-keeping', the way this is used in the Jacobi reduction of the Hamiltonian, the chain of canonical transformations leading the Hamiltonian to its final form in Poincaré canonical variables, as well as a number of precision tests about the order of the multipole truncation, comparison with the Laplace-Lagrange series, etc.  Section <ref> analyzes the dynamics and the phase portraits under the secular Hamiltonian model computed as in section <ref>. In particular, subsections <ref> and <ref> discuss the phase-space properties of the integrable approximations $\Hscr_{planar}$ and $\Hscr_{int}$, and introduces the basic terminology for the orbits of the apsidal corotation type (modes A and B). Subsection <ref> deals, instead,  with an analysis of the phase-space structure under the full (non-integrable) Hamiltonian (<ref>).  This includes the semi-analytical method (via a normal form) used to demonstrate the continuation of the A and B modes as well as to produce the time series yielding the evolution of the eccentricity vectors for both planets along these modes. Subsection <ref> discusses numerical results obtained by computing the surfaces of section as the energy $\Escr$ increases (or the parameter $\delta \Escr$ decreases), showing the sequence of bifurcations that connect the apsidal corotation with the Lidov-Kozai regime. Subsection <ref> contains some results related to the `Kozai mechanism', i.e.,  the transition from linear stability to instability for the inclined circular orbit of one of the two planets. These results are compared to an analytical approximation obtained in the framework of the quadrupolar approximation. Section <ref> contains the investigation of how the phase portraits change in cases with different distance or mass ratios of the planets, covering various hierachical models ($r_2/r_3=1/7,\,1/3,\,$ $m_2/m_3=1/10,\,1/3,\,1,\,3,\,10\,$). Section <ref> provides a summary of the main conclusions from the present study.

§ HAMILTONIAN MODEL

§.§ Averaged Hamiltonian

We will focus below mostly on the properties of a secular model $\Hscr_{sec}$ for the Hamiltonian (<ref>) obtained by performing averaging with respect to the fast angles just `by scissors'. We denote by $(a,\e,\i,M,\omega,\Omega)$ the Keplerian elements of a body (semi-major axis, eccentricity, inclination, mean anomaly, argument of the periastron, argument of the nodes), and by $\lambda=M+\omega+\Omega$, $\varpi=\omega+\Omega$ the mean longitude and longitude of the periastron respectively. We then have the following properties of the secular model obtained after averaging `by scissors':

\begin{eqnarray}\label{hamscis}
\Hscr_{sec}&=&{1\over 4\pi^2}
\int_{0}^{2\pi}\int_{0}^{2\pi} H(\vet{r_2},\vet{r_3},\vet{p_2},\vet{p_3}) \,dM_{2}\,dM_{3} \nonumber\\
{1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}
\left(\left({1\over m_0}+{1\over m_2}\right)\frac{\vet{p_2}^2}{2}-\frac{\mathcal{G}\, m_0\, m_2}{r_2}
\left({1\over m_0}+{1\over m_3}\right)\frac{\vet{p_3}^2}{2}-\frac{\mathcal{G}\, m_0\, m_3}{r_3}
\right)\,dM_{2}\,dM_{3} \nonumber\\
&+&{1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}
\frac{\vet{p_2}\cdot\vet{p_3}}{ m_0}dM_{2}\,dM_{3}
-{1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{\mathcal{G}\, m_2\, m_3}{|\vet{r_2}-\vet{r_3}|}\,dM_{2}\,dM_{3}~~.
\end{eqnarray}

We have

\begin{eqnarray}\label{hkep0sec}
~&~& {1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}
\left(\left({1\over m_0}+{1\over m_2}\right)\frac{\vet{p_2}^2}{2}-\frac{\mathcal{G}\, m_0\, m_2}{r_2}
\left({1\over m_0}+{1\over m_3}\right)\frac{\vet{p_3}^2}{2}-\frac{\mathcal{G}\, m_0\, m_3}{r_3}
\right)\,dM_{2}\,dM_{3}  \nonumber\\
&=&-{\mathcal{G} m_0 m_2\over 2a_2}-{\mathcal{G} m_0 m_3\over 2a_3}+
{\mathcal{G}  m_2^2\over 2a_2}+{\mathcal{G} m_3^2\over 2a_3}~~.
\end{eqnarray}

The indirect part of the disturbing function, depending on the product $\vet{p_2}\cdot\vet{p_3}$, yields
zero average

\begin{equation}\label{indmomenta}
{1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{\vet{p_2}\cdot\vet{p_3}}{ m_0}\,dM_{2}\,dM_{3}
={1\over 4\pi^2}
{\Gscr \,m_2 m_3\over \sqrt{a_2^3\,a_3^3}}
\left( \int_{0}^{2\pi}{d\vet{r_2}\over dM_2}  \,dM_{2}
\int_{0}^{2\pi}{d\vet{r_3}\over dM_3}  \,dM_{3}\right)=0~~.
\end{equation}

To compute the average of the direct part $\mathcal{G}m_2m_3/|\vet{r_2}-\vet{r_3}|$, we assume a dynamical regime of the planetary system in which the distance ratio $r_2/r_3$ (where $r_{i}=|\vet{r}_{i}|\,$, $i=2,3$) remains always smaller than unity. Then, the quantity $\mathcal{G}m_2m_3/|\vet{r_2}-\vet{r_3}|$ admits a convergent Legendre multipolar expansion in powers of the quantity $r_2/r_3<1$:

\begin{align*}
-{1\over 4\pi^2}\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{\mathcal{G}\, m_2\, m_3}{|\vet{r_2}-\vet{r_3}|} \,&dM_{2}\,dM_{3}\, =\\
&-{\mathcal{G}\, m_2\, m_3\over 4\pi^2}
\int_{0}^{2\pi}\int_{0}^{2\pi}\left(\frac{1}{r_{3}} +\frac{\vet{r}_{2}\cdot\vet{r}_{3}}{r_{3}^3}-\frac{1}{2}\frac{r_2^2}{r_3^3}+\frac{3}{2}\frac{(\vet{r}_{2}\cdot\vet{r}_{3})^2}{r_{3}^5} +\ldots\right)
\,dM_{2}\,dM_{3}\, .
\end{align*}

We have:

\begin{align*}
&-{\mathcal{G}\, m_2\, m_3\over 4\pi^2}
\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{1}{r_{3}} dM_1dM_2=-\frac{\mathcal{G}m_2 m_3}{a_3}\, , &
\end{align*}

We then need to compute

\begin{equation}\label{remainder}
\mathcal{R}_{sec}= {1\over 4\pi^2}
\int_{0}^{2\pi}\int_{0}^{2\pi}-\frac{\mathcal{G}\, m_2\, m_3}{r_3}\left(-\frac{1}{2}\frac{r_2^2}{r_3^2}+\frac{3}{2}\frac{(\vet{r}_{2}\cdot\vet{r}_{3})^2}{r_{3}^4} +\ldots\right) \,dM_{2}\,dM_{3}\, .
\end{equation}

Keeping only the lowest order term in the integrand of (<ref>), proportional to $(r_2/r_3)^2\,$, leads to the so-called quadrupole approximation; the next truncation (up to terms proportional to $(r_2/r_3)^3\,$) is the octupole approximation, etc. The integrals of any multipole approximation  can be computed in so-called closed form (i.e. without expansions in the eccentricities), by avoiding completely the series reversion of Kepler's equation, using, instead, the change of variables $M_2\rightarrow u_2$ (eccentric anomaly), $M_3\rightarrow f_3\,$ (true anomaly). We have

\begin{equation}\label{trick.closed}
dM_2=(1-\e_2\cos u_2)\,du_2,~dM_3={r_3^2\over a_3^2\sqrt{1-\e_3^2}} \, df_3,~
r_2=a_2(1-\e_2\cos(u_2)),~\frac{1}{r_3}=\frac{1+\e_3\cos(f_3)}{a_3\,(1-\e_3^2)}\,  .
\end{equation}

Replacing the above expressions in (<ref>) and performing all trigonometric reductions, we find a trigonometric polynomial series containing only terms of the form $\cos(K_2u_2+K_3f_3+...)$, with $K_2,K_3$ integers. This implies that the average can be computed by just scissor-cutting all the terms in the integrand of  (<ref>) for which $|K_2|+|K_3|\neq 0$. This leads to a closed-form expression for the secular Hamiltonian

\begin{equation}\label{H.closed}
\Hscr_{sec}=-{\mathcal{G} m_0 m_2\over 2a_2}-{\mathcal{G} m_0 m_3\over 2a_3}+{\mathcal{G}  m_2^2\over 2a_2}+{\mathcal{G} m_3^2\over 2a_3}-{\mathcal{G} m_2 m_3\over a_3}+\mathcal{R}_{sec}(a_2,a_3,\e_2,\e_3,\i_2,\i_3,\omega_2,\omega_3,\Omega_2-\Omega_3)~~.
\end{equation}

The following are some relevant remarks regarding the Hamiltonian (<ref>):

Remark 1: the averaging (<ref>) yields a valid secular model only when the system is assumed to be far from any low-order mean-motion resonance. By `low' it is implied that no resonance condition of the form

should be satisfied, with $K_2,K_3$ integers, $\mu=\max(m_2/m_0,m_3/m_0)$, $n_2=\left(\mathcal{G}(m_0+m_2)/a_2^3\right)^{1/2}$, $n_3=\left(\mathcal{G}(m_0+m_3)/a_3^3\right)^{1/2}$, for an order $|K_2|+|K_3|$ inferior or equal to the order $N_P$ of the multipole expansion in Eq.(<ref>). Formally, such a requirement reflects the fact that the averaging by scissors serves to substitute the complete procedure of first-order averaging, which involves a canonical transformation to properly eliminate from the Hamiltonian the fast angles $M_2,M_3$. Such a transformation can be defined in closed form (see, for example, [Cavallari and Efthymiopoulos, 2022]), but it involves divisors of the form $K_2n_2+K_3n_3$ which become very small near any mean motion resonance of order smaller or equal to $N_P$.

Remark 2: the function $\mathcal{R}_{sec}$ is trigonometric polynomial in the angles $\omega_2,\omega_3$, as well as $\Omega_2-\Omega_3$. The dependence on $\Omega_2,\Omega_3$ only by the difference $\Omega_2-\Omega_3$ is a consequence of the fact that the sum of the angular momenta $H_2+H_3$ perpendicularly to the system's Laplace plane is an invariant of the motion of the complete Hamiltonian (<ref>), which is preserved also in the secular Hamiltonian (<ref>).

Remark 3: the function $\mathcal{R}_{sec}$ is polynomial in the quantities $a_2,\e_2,\e_3,\sin \i_2,\sin \i_3,\cos \i_2,\cos \i_3$, while it is rational in the quantities $a_3,\eta_3=\sqrt{1-\e_3^2}$. As regards the dependence on $\sin \i_2,\sin \i_3,\cos \i_2,\cos \i_3$, this follows the symmetry that $\mathcal{R}_{sec}$ can only depend on the quantity $\cos(\i_2+\i_3)$, where $\i_2+\i_3=\i_{mut}$ is the mutual inclination of the planetary trajectories. This symmetry is exploited when performing Jacobi reduction of the nodes in the Hamiltonian $\Hscr_{sec}$ (see next subsection). On the other hand, the rational dependence on the quantity $\sqrt{1-\e_3^2}$ implies that any expansion in the eccentricities performed after the closed-form averaging is convergent in the entire domain $|\e_2|<1,|\e_3|<1$. This allowance to perform an a posteriori expansion of the Hamiltonian $\Hscr_{sec}$ in the planetary eccentricities proves useful in the theoretical analysis of the dynamical properties of such a Hamiltonian (see section <ref>).

§.§ Book-keeping and Jacobi reduction of the nodes

Consider the set of Delaunay canonical variables:

\begin{align}\label{Delau}
&L_j=m_j\sqrt{\mathcal{G}\,m_0\,a_j}\, , & &l_j=M_j\, ,\notag\\
&G_j=L_j\sqrt{1-\e^2_j}\, , & &g_j=\omega_j\, , \\
&H_j=G_j\cos(\i_j)\, , & &h_j=\Omega_j\, \notag
\end{align}

with $j=2,3$. Since $a_j=a_j(L_j)$, $\e_j=\e_j(L_j,G_j)$, $\i_j=\i_j(L_j,G_j,H_j)$, substituting the corresponding expressions in the Hamiltonian $\mathcal{H}_{sec}$ (Eq. (<ref>)) leads to a Hamiltonian model with six degrees of freedom

\begin{equation}\label{hsecdelau}
\mathcal{H}_{\sec}\equiv\mathcal{H}_{sec}(L_2,L_3,G_2,G_3,H_2,H_3,\omega_2,\omega_3,\Omega_2-\Omega_3)~~.
\end{equation}

Since the angles $M_2,M_3$ are ignorable, the Delaunay momenta $L_2,L_3$ are integrals of motion under the Hamiltonian flow of $\mathcal{H}_{sec}$. There are two more integrals of motion in involution with $L_2,L_3$, corresponding to the two components of the total angular momentum, normal and parallel to the Laplace plane. For the normal component we have $L_z=H_2+H_3=C=const$, while for the parallel component we have $L_{\parallel}=G_2\sin(\i_2(L_2,G_2,H_2))-G_3\sin(\i_3(L_3,G_3,H_3))=0$. Thus, with an appropriate reduction, we can reduce the system to a Hamiltonian with two degrees of freedom. This reduction, implemented on the invariant manifold $L_z=C$, $L_\parallel=0$ is called the Jacobi reduction of the nodes [Jacobi, 1842] and it exploits one more invariance of the system: since $L_\parallel=0$ we necessarily have $\Omega_2-\Omega_3=\pi=const$.

In our analytical treatment of the Hamiltonian $\mathcal{H}_{sec}$ it turns convenient to perform the Jacobi reduction of the nodes by simultaneously introducing two `book-keeping symbols', $\varepsilon$, $\eta$ [Efthymiopoulos, 2012], both with numerical values $\varepsilon=1$, $\eta=1$, whose role is the following: i) $\varepsilon$ keeps track of the order of a certain term in the planetary eccentricities and inclinations, ii) $\eta$ separates Hamiltonian terms which depend on powers of the quantity $\cos(\i_2+\i_3)$ from those terms which do not depend on the mutual inclination $\i_{mut}=\i_2+\i_3$.

More specifically, we Jacobi-reduce the Hamiltonian $\mathcal{H}_{sec}$ by the following steps:

Step 1: Canonical transformation. Similarly as in [Libert and Henrard, 2007], we introduce the canonical transformation

\begin{align}\label{Delaujac}
&\Lambda_2=L_2\, , & &\lambda_2=M_2+\omega_2+\Omega_2\, ,\notag\\
&\Lambda_3=L_3\, , & &\lambda_3=M_3+\omega_3+\Omega_3\, ,\notag\\
&W_2=L_2-G_2\,   , & &w_2=-\omega_2\, , \\
&W_3=L_3-G_3\,   , & &w_3=-\omega_3\, , \notag\\
&R_2=L_2-H_2\,   , & &\theta_{r2}=\Omega_3-\Omega_2\, , \notag\\
&R_3=L_2+L_3-H_2-H_3=\AMD\,   , & &\theta_{r3}=-\Omega_3\, . \notag
\end{align}

Inverting the transformation (<ref>) and substituting the result into the Hamiltonian (<ref>), the secular Hamiltonian obtains the form:

\begin{equation}\label{hsecdelaujac}
\mathcal{H}_{\sec}\equiv
\mathcal{H}_{sec}(\Lambda_2,\Lambda_3,W_2,W_3,R_2,R_3,w_2,w_3,\theta_{r2})~~.
\end{equation}

Step 2: Jacobi reduction with book-keeping. Hamilton's equations for the Hamiltonian (<ref>) yield $\dot{\Lambda}_2=\dot{\Lambda}_3=0$ (constancy of the semi-major axes $a_2$, $a_3$), as well as $\dot{R}_3=0$ (constancy of $L_z$, implying the constancy of the AMD). The crucial point stems from the following invariance property of
the Hamiltonian (<ref>):

\begin{equation}\label{dhamdthr2}
\dot{\theta}_{r2}=
{\partial\mathcal{H}_{sec}(\Lambda_2,\Lambda_3,W_2,W_3,R_2,R_3=\AMD,w_2,w_3,\theta_{r2}=\pi)\over\partial R_2}=0
\end{equation}

corresponding to the invariance in time of the relation $\theta_{r2}=\Omega_3-\Omega_2=\pi$ for all trajectories. Equation (<ref>), however, implies that when the substitution $\Omega_3=\Omega_2+\pi$ is made in the Hamiltonian (<ref>), the resulting expression becomes independent of $R_2$, and the reduced set of Hamilton's equations

\begin{align}\label{hameqred}
\over\partial W_2}
\, , \notag\\
\over\partial W_3}
\, , \\
\over\partial w_2}
\, , \notag\\
\over\partial w_3}
\, , \notag
\end{align}

remains valid. Note that $\Lambda_2,\,\Lambda_3$ in the above secular equations of motion are constant, and can be effectively treated as parameters, depending on the (also constant) parameters $a_2,a_3$. Also, the terms

-{\mathcal{G} m_0 m_2\over 2a_2}-{\mathcal{G} m_0 m_3\over 2a_3}+{\mathcal{G}  m_2^2\over 2a_2}+{\mathcal{G} m_3^2\over 2a_3}-{\mathcal{G} m_2 m_3\over a_3}

appearing in Eq. (<ref>), do not contribute to the secular equations of motion, and can be omitted from further analysis, by just renaming $\mathcal{R}_{sec}\rightarrow\mathcal{H}_{sec}$. Finally, Eqs. (<ref>) allow to compute the secular evolution of only the eccentricity vectors of the two planets. To obtain the evolution in inclination, instead, we use the following relations,  obtained directly from the conservation of the angular momentum:

\begin{equation}\label{cosi23e23}
\cos(\i_2)={\Lambda_2^2(1-\e_2^2)-\Lambda_3^2(1-\e_3^2)+L_z^2\over 2L_z\Lambda_2\sqrt{1-\e_2^2}},~~~
\cos(\i_3)={L_z-\Lambda_2\sqrt{1-\e_2^2}\cos(\i_2)\over \Lambda_3\sqrt{1-\e_3^2}}~~.
\end{equation}

In the practical implementation of step 2, we perform the substitutions of the angle variables $\Omega_3=\Omega_2+\pi$, as well as $\omega_2=-w_2$, $\omega_3=-w_3$ in the Hamiltonian, but leave implicit the latter's dependence on the action variables $(\Lambda_2,\Lambda_3,W_2,W_2,R_2,R_3)$ through the elements $(a_2,a_3,\e_2,\e_3,\i_2,\i_3)$. Owing, however, to the non-dependence of the Hamiltonian on $R_2$ after the substitution $\Omega_3=\Omega_2+\pi$, we have that the Hamiltonian depends on the inclinations only through the trigonometric combination $\cos(\i_2+\i_3)=\cos(\i_2)\cos(\i_3)-\sin(\i_2)\sin(\i_3)\,$. Taking into account that

\left(1- \cos(\i_2)\cos(\i_3)\right)=\Oscr(\varepsilon^2)\, , \qquad\,
\sin(\i_2)\sin(\i_3)=\Oscr(\varepsilon^2)~~,

where the attribution of $\Oscr(\varepsilon^2)$ to an expression stands for `second order in the eccentricities and inclinations', we then introduce the following book-keeping control identities:

\begin{align}\label{riduzione.cos.sen}
&\cos(\i_2)\cos(\i_3)=\varepsilon^2\eta\left(\cos(\i_2)\cos(\i_3)-1\right)+1\, , &
&\sin(\i_2)\sin(\i_3)=\varepsilon^2\eta\sin(\i_2)\sin(\i_3)\, .
\end{align}

We also use the substitution rule

\begin{equation}\label{seni2i3}
\sin(\i_2)\sin(\i_3)= \cos(\i_2)\cos(\i_3)+\frac{\Lambda_2\sqrt{1-\e_2^2}}{2\,\Lambda_3\sqrt{1-\e_3^2}}
-\frac{L_z^2}{2\,\Lambda_2\sqrt{1-\e_2^2}\,\Lambda_3\sqrt{1-\e_3^2}}\, .
\end{equation}

Substituting the above expressions into the Hamiltonian, and truncating the resulting expression up to a pre-selected maximum order in book-keeping $N_{bk}$, by symmetry the terms in equal powers of the products $\sin(\i_2)\sin(\i_3)$ and $\cos(\i_2)\cos(\i_3)$ are opposite, and thus they are cancelled. Hence, the Hamiltonian resumes the form:

\begin{equation}\label{hambketa}
\mathcal{H}_{sec}=\sum_{s_1=0}^{N_{bk}/2}\eta^{s_1}\varepsilon^{2s_1}\mathcal{H}_{sec,s1}(\e_2,\e_3,w_2,w_3;a_2,a_3,L_z)~~.
\end{equation}

We finally restore the values $\varepsilon=1$ and $\eta=1$ of the book-keeping parameters, and write the truncated (up to book-keeping order $N_{bk}$) Hamiltonian as:

\begin{equation}\label{hamplsp}
\mathcal{H}_{sec}=\mathcal{H}_{planar}+\mathcal{H}_{space}
\end{equation}

\begin{equation}\label{hamplspanal}
\mathcal{H}_{planar}= \mathcal{H}_{sec,0}(\e_2,\e_3,w_2-w_3;a_2,a_3),~~~
\mathcal{H}_{space}= \sum_{s_1=1}^{N_{bk}/2}\mathcal{H}_{sec,s1}(\e_2,\e_3,w_2,w_3;a_2,a_3,L_z)~.
\end{equation}

Note that the term $\mathcal{H}_{planar}$ depends only on the difference $w_2-w_3$, since, in the planar case, the sum $W_2+W_3$, which represents the total angular momentum deficit for planar orbits, is a conserved quantity.

Step 3: Expansion in the orbital eccentricities. We re-introduce the book-keeping:

\begin{equation}\label{ecclzbk}
\e_2\rightarrow \varepsilon \e_2,~~\e_3\rightarrow \varepsilon \e_3,~~
L_z \rightarrow \Lambda_2+\Lambda_3-\varepsilon^2\AMD\, .
\end{equation}

Substituting the above expressions in the Hamiltonian (<ref>), and expanding the resulting expressions in powers of the book-keeping parameter $\varepsilon$, apart from a term depending only on the constants $a_1,a_2$ and $\AMD$, we arrive at the truncated Hamiltonian:

\begin{eqnarray}\label{hamnostra}
\label{Ham.Lapl.nostra2}
\Hscr_{sec}&=&\Hscr_{planar}+\Hscr_{space} \\
&=&\sum_{s_2=1}^{N_{bk}/2}\varepsilon^{2 s_{2}} h_{s_{2}}(\e_2,\e_3, w_2-w_3;a_2,a_3)
+\sum_{s_2=1}^{N_{bk}/2}\varepsilon^{2 s_{2}}\tilde{h}_{s_2}(\e_2,\e_3,w_2,w_3;a_2,a_3,\AMD)~~.\nonumber
\end{eqnarray}

Step 4: introduction of Poincaré variables:
The Hamiltonian (<ref>) is polynomial in the eccentricities $\e_2,\,\e_3$. It was already stressed that the series expansion with respect to the quantities $\eta_2=\sqrt{1-\e_2^2}$ and $\eta_3=\sqrt{1-\e_3^2}$ introduces no divergence in the polydisc $|\e_2|<1$, $|\e_3|<1$. Besides, the Hamiltonian satisfies the D'Alembert rules in the eccentricities: every trigonometric term in it is of the form $h_{\ell_2,\ell_3,k_2,k_3}(a_2,a_3)(\AMD)^{\ell_1}\e_2^{\ell_2}\e_3^{\ell_3}{\cos}(k_2w_2+k_3w_3)$, with $\ell_1$ positive integer, and the integers $\ell_2,\ell_3,k_2,k_3$ satisfying i) $\ell_2,\ell_3>0$, ii) $\ell_2+\ell_3\geq |k_2|+|k_3|$, and iii) $mod(\ell_2,2)=mod(k_2,2)$, $mod(\ell_3,2)=mod(k_3,2)$. Rules (i) to (iii) imply, now, that the Hamiltonian is polynomial in the Poincaré canonical variables stemming from the variables $(w_j,W_j)$, $j=2,3$ through the canonical transformation

\begin{align}\label{Poincare1}
& X_j= -\sqrt{2 \,W_j}\cos(w_j)\, , & & Y_j=\sqrt{2\,W_j}\sin(w_j)\, , & &j=2,3~~.
\end{align}

To obtain a truncated polynomial series of the Hamiltonian in the variables $(X_j,Y_j)$, we substitute into the Hamiltonian (<ref>) the expressions

\sin(w_j)={Y_j\over\sqrt{2W_j}},~~~\cos(w_j)=-{X_j\over\sqrt{2W_j}},~~~\e_j=\sqrt{{2W_j\over\Lambda_j}-\varepsilon^2\left({W_j\over\Lambda_j}\right)^2},~~j=2,3~

and expand the result in powers of the book-keeping parameter $\varepsilon$, up to the truncation order $N_{bk}$. This leads to an expression which no longer contains trigonometric functions of the angles $w_j$, while it still contains integer powers of the actions $W_j$. We then substitute $W_j\rightarrow(X_j^2+Y_j^2)/2$, and, finally, set back $\varepsilon=1$. In this way we arrive at the final secular model:

\begin{eqnarray}\label{Ham.final}
\Hscr_{sec}(X_2, X_3, Y_2, Y_3; \AMD)&=&
\Hscr_{planar}(X_2, X_3, Y_2, Y_3)+\Hscr_{space}(X_2, X_3, Y_2, Y_3;\AMD) \nonumber\\
\sum_{\ell\in\mathbb{N}^4,|\ell|=2}^{N_{bk}}
\mathcal{K}_{plane,\ell}X_2^{\ell_1}\,X_3^{\ell_2}\,Y_2^{\ell_3}\,Y_3^{\ell_4} \\
\sum_{\ell\in\mathbb{N}^4,|\ell|=2}^{N_{bk}}
\mathcal{K}_{space,\ell}(\AMD)X_2^{\ell_1}\,X_3^{\ell_2}\,Y_2^{\ell_3}\,Y_3^{\ell_4}~~.\nonumber
\end{eqnarray}

Note that, by symmetry, the value of the disturbing function (and hence of the secular Hamiltonian) remains invariant by the rotation of both planets' argument of the pericenter by $\pi$, hence the secular Hamiltonian is necessarily even in the planetary eccentricities. This implies that the Hamiltonian (<ref>) contains only even powers of the Poincaré variables $(X_j,Y_j)\,$, i.e. $\ell_1+\ell_2$ and $\ell_3+\ell_4$ are even. In particular, Hamilton's equations take the form:

\begin{equation}\label{conseguenza.Lib-Henr}
\begin{split}
&\dot{X_i}=\frac{\partial \Hscr_{sec}}{\partial Y_i}
=Y_2\,F_{X_i,2}(X_2, X_3, Y_2, Y_3)+Y_3\,F_{X_i,3}(X_2, X_3, Y_2, Y_3)\, ,\\
&\dot{Y_i}=-\frac{\partial \Hscr_{sec}}{\partial X_i}
=X_2\,F_{Y_i,2}(X_2, X_3, Y_2, Y_3)+X_3\,F_{Y_i,3}(X_2, X_3, Y_2, Y_3)\, ,
\end{split}
\end{equation}

where the polynomials $F_{X_i,2}$, $F_{X_i,3}$, $F_{Y_i,2}$, $F_{Y_i,2}\,$, $ i=2,\,3\, ,$ are even, starting with constant terms.

§.§ Poincaré surface of section. Precision tests

We will now discuss several precision tests, based on the method of comparison of phase portraits, which aim to establish which is the minimum multipole order, as well as minimum order in the eccentricities of the Hamiltonian $\Hscr_{sec}$, such that the Hamiltonian, truncated at the above orders, represents with sufficient precision the dynamics at the timescales of the secular system produced by averaging of the original Hamiltonian (<ref>).

§.§.§ Poincaré surface of section: definitions

To visualize phase portraits, use is made below of a Poincaré surface of section $\mathcal{P}(\Escr;\AMD)$ defined by the relations:

\begin{equation}\label{pcsec}
\begin{aligned}
\mathcal{P}(\Escr;\AMD)=&
\bigg\{ (X_2,Y_2,X_3,Y_3)\in\mathbb{R}^4:~
Y_3=0, ~
\dot{Y}_3=-{\partial\Hscr_{sec}(X_2,Y_2,X_3,Y_3=0;\AMD)\over\partial X_3}\geq 0,~\\
\cos(\i_{max})\leq\cos(\i_{mut})(X_2,Y_2,X_3,Y_3=0;\AMD)\leq 1\bigg\}~~,
\end{aligned}
\end{equation}

where the mutual inclination $\i_{mut}=\i_2+\i_3$ for fixed $\AMD$, or angular momentum $L_z=\Lambda_2+\Lambda_3-\AMD\,$, is given by

\begin{equation}\label{cosimutua}
\cos(\i_{mut})=
\frac{L_z^2-\Lambda_2^2-\Lambda_3^2+\Lambda_2^2\,\e_2^2+\Lambda_3^2\,\e_3^2}{2\Lambda_2\Lambda_3\sqrt{1-\e_2^2}\sqrt{1-\e_3^2}}
\frac{L_z^2-\Lambda_2^2\left(1-{X_2^2+Y_2^2\over 2L_2}\right)^2
-\Lambda_3^2\left(1-{X_3^2+Y_3^2\over 2L_3}\right)^2}
{2\Lambda_2\Lambda_3\left(1-{X_2^2+Y_2^2\over 2L_2}\right)\left(1-{X_3^2+Y_3^2\over 2L_3}\right)}
\end{equation}

and the maximum possible mutual inclination consistent with the given $\AMD$ is

\begin{equation}\label{imutuamax}
\i_{max}=
\cos^{-1}\left(\frac{L_z^2-\Lambda_2^2-\Lambda_3^2}{2\Lambda_2\Lambda_3}\right)~~.
\end{equation}

The phase portrait corresponding to the Poincaré surface of section at a fixed level of energy $\Escr$ is obtained numerically, by choosing several initial conditions $(X_2,Y_2)\in\mathcal{D}(\Escr)\subset\mathbb{R}^2$, where $\mathcal{D}(\Escr)$ is the domain of permissible initial conditions consistent with the definition of the surface of section as in Eq.(<ref>). For each initial condition, we then iterate the corresponding orbit under the Hamiltonian $\Hscr_{sec}$, and plot the consequent points $(X_2,Y_2)$, or, equivalently $(\e_2\cos\omega_2,\e_2\sin\omega_2)$, with

\begin{equation}\label{xy2e2ome}
\e_2=\left(1-(1-(X_2^2+Y_2^2)/(2\Lambda_2))^2\right)^{1/2},~~~\omega_2=\mathrm{sgn}\left(\frac{-Y_2}{\sqrt{X_2^2+Y_2^2}}\right)\arccos\left(\frac{-X_2}{\sqrt{X_2^2+Y_2^2}}\right)
\end{equation}

every time when the orbit intersects the surface of section. Note that the conditions $Y_3=0,\dot{Y}_3\geq 0$ imply i) that the Poincaré mapping defined by the successive iterates is symplectic, and ii) that the section definition corresponds physically to instants when the orbit of the outer planet crosses the pericenter $\omega_3=\pi$. The symplecticity of the Poincaré mapping, along with the straightforward physical interpretation, are the main motives for the choice of variables and for the visualization of phase portraits via the definition of the Poincaré section as above.

The domain $\mathcal{D}(\Escr)$ is non-null in a range of energies $\Escr_{min}\leq\Escr<0$. The energy $\Escr_{min}$ is computed as follows: consider the surface $\mathcal{I}_0$ of all possible points satisfying the section condition $Y_3=0$ as well as the lowermost limit of possible mutual inclination $\i_{mut}=0$. The surface  $\mathcal{I}_0$ is the sphere given by

\begin{equation}\label{imut0}
\mathcal{I}_0=\bigg\{X_3^2+X_2^2+Y_2^2=2\AMD\bigg\}~~.
\end{equation}

We then find the point $(X_{2,0},Y_{2,0})$ on the surface $\mathcal{I}_0$ where the energy is minimum by looking to the solutions of the system of equations

\begin{eqnarray}\label{eqenmin}
{\partial\Hscr_{sec}(X_2,Y_2,X_3^2=2\AMD-X_2^2-Y_2^2,Y_3=0; \AMD)\over\partial X_2}=0 \\
{\partial\Hscr_{sec}(X_2,Y_2,X_3^2=2\AMD-X_2^2-Y_2^2,Y_3=0; \AMD)\over\partial Y_2}=0~,\nonumber
\end{eqnarray}

satisfying also the section condition

\begin{equation}\label{y3dotsec}
\dot{Y}_3=\left({\partial\Hscr_{sec}(X_2,Y_2,X_3,Y_3=0;\AMD)\over\partial X_3}\right)_{X_3^2=2\AMD-X_2^2-Y_2^2}\geq 0~~.
\end{equation}

We find two solutions $(X_{2,0},Y_{2,0}=0)\,$ and $(X_{2,1},Y_{2,1}=0)$ of the system of equations (<ref>) and (<ref>); however, only $(X_{2,0},Y_{2,0}=0)\,$ corresponds to the minimum value of the energy, i.e. $\Hscr(X_2=X_{2,0}, Y_2=Y_{2,0}=0, X_3= X_{3,0}=\sqrt{2\AMD - X_{2,0}^2}, Y_3=0 )=\Escr_{min}$, as verified by the Hessian matrix

\begin{pmatrix}
&\displaystyle{\frac{\partial^2 \Hscr_{sec}(X_3=\sqrt{2\AMD - X_{2}^2},Y_3=0)}{\partial X_2^2}}
& \displaystyle{\frac{\partial^2 \Hscr_{sec}(X_3=\sqrt{2\AMD - X_{2}^2},Y_3=0)}{\partial Y_2\partial X_2}}
\\
&\displaystyle{\frac{\partial^2 \Hscr_{sec}(X_3=\sqrt{2\AMD - X_{2}^2},Y_3=0)}{\partial X_2\partial Y_2}}
& \displaystyle{\frac{\partial^2 \Hscr_{sec}(X_3=\sqrt{2\AMD - X_{2}^2},Y_3=0)}{\partial Y_2^2}}
\end{pmatrix}_{(X_2=X_{2,0}, Y_2=Y_{2,0})}

which is positive definite. On the other hand, the solution $(X_{2,1},\,Y_{2,1}=0,\,X_{3,1}=\sqrt{2\AMD - X_{2,1}^2})$ corresponds to the energy $\Escr_{2,3}$, and yield a negative definite Hessian matrix. Thus, $\Escr_{2,3}$ corresponds to the maximum energy for which we can have co-planar orbits (see discussion below).

The uniqueness of the solution realizing the minimal value of the energy can be understood by the following argument: for every fixed value of the energy $\Escr$, sufficiently close to the origin the restriction of the manifold of constant energy to the section $Y_3=0$

\begin{equation}\label{manifene}
\mathcal{M}(\Escr)=\bigg\{(X_2,Y_2,X_3)\in\mathbb{R}^3: \Hscr_{sec}(X_2,X_3,Y_2,Y_3=0;\AMD)=\Escr\bigg\}
\end{equation}

yields a surface. We will show in the next section that, for energies smaller than a suitably defined threshold,  the surface $\mathcal{M}(\Escr)$ is close to the surface of an integrable model

\begin{equation}\label{manifeneint}
\mathcal{M}_{int}(\Escr)=\bigg\{(X_2,Y_2,X_3)\in\mathbb{R}^3: \Hscr_{int}(X_2,X_3,Y_2,Y_3=0;\AMD)=\Escr\bigg\}
\end{equation}

where $\Hscr_{int}$ contains $\Hscr_{planar}$ as well as a part of $\Hscr_{space}$. The manifolds $\mathcal{M}_{int}(\Escr)$, in turn, are ellipsoidal-like closed convex surfaces. The ellipsoidal form follows from the fact that the leading order terms of the restriction of $\Hscr_{int}$ to the section $Y_3=0$ yield a negative-definite quadratic form

\Hscr_{int}(X_2,X_3,Y_2,Y_3=0)=-\nu_2\left({X_2^2+Y_2^2\over 2}\right)-\nu_3{X_3^2\over 2}+c_{23}X_2X_3+\ldots

with $\nu_2,\nu_3>0$, $|c_{2,3}|<\min(\nu_2,\nu_3)$.

Owing to their proximity to the manifolds $\mathcal{M}_{int}(\Escr)$, for energies low enough, the manifolds $\mathcal{M}(\Escr)$ are also ellipsoidal-like (the exact evolution of the form of $\mathcal{M}(\Escr)$ as the energy increases will be discussed in the next section). As a consequence, there is a negative value of the energy $\Escr=\Escr_{min}$ such that for $\Escr<\Escr_{min}$ the ellipsoidal manifold $\mathcal{M}(\Escr)$ surrounds the sphere $\mathcal{I}_0$, thus leading to the unphysical condition $\cos(\i_{mut})>1$. At the energy $\Escr=\Escr_{min}$ the ellipsoidal manifold $\mathcal{M}(\Escr_{min})$ has two points of tangency with the sphere $\mathcal{I}_0$ (see Fig.<ref>). Both tangencies occur at the plane $Y_2=0$, but only one of them satisfies the condition $\dot{Y}_3\geq 0$. Then, we set $(X_{2,0},Y_{2,0}=0,X_{3,0}^2=2\AMD-X_{2,0}^2\,)$ equal to the coordinates of the corresponding point of tangency, and $\Escr_{min}=\Hscr_{sec}(X_2=X_{2,0},X_3=X_{3,0},Y_2=0,Y_3=0;\AMD)$.

Left: tangency between the ellipsoidal surfaces of minimum possible energy $\mathcal{M}(\Escr_{min})$, with $\Escr_{min}=-1.18237\cdot 10^{-4}$ (light gray) and the sphere $\mathcal{I}_0$ (gray). The only permissible initial condition on the surface of section corresponds to the point of tangency (black dot). The green slightly inclined surface  represents the condition $\dot{Y}_3=0$. Right: the tangency between the surfaces $\mathcal{M}(\Escr_{min})$  (thick light gray closed curves) and the sphere $\mathcal{I}_0$ (thin gray curve) as seen in the plane $(X_2,X_3)$ for $Y_2=0$. The green slightly inclined curve represents the condition $\dot{Y}_3=0$ and separates the plane $(X_2,X_3)$ in an upper domain, where $\dot{Y}_3>0$, and a lower domain ($\dot{Y}_3<0$). The point of tangency (thick dot) marks a fixed point corresponding to a planar periodic orbit of the apsidal corotation type called `mode A' (anti-aligned) in section <ref>, yielding $(X_2 = X_{2,0}=0.101237$, $X_3 = X_{3,0}=\sqrt{2\AMD-X_{2,0}^2}=0.148862$, $Y_2 = Y_{2,0}=0$, $Y_3=0\,)$, that correspond to $\e_2=0.337415$, $\e_3=0.433811$.

In order to numerically specify, now, the limits of the domain $\mathcal{D}(\Escr)$ on the Poincaré section for any value of the energy in the range $\Escr_{min}<\Escr<0$ we work as follows: fixing any value of the angle $w_2=-\omega_2$ in the interval $0\leq w_2\leq\pi$, the line defined parametrically by the relations:

\begin{equation}\label{linexy2s2}
\mathcal{L}(w_2):\{X_2=s_2\cos(w_2),~Y_2=s_2\sin(w_2),s_2\in\mathbb{R}\}
\end{equation}

defines a plane $\mathcal{PL}(w_2):\{(X_2,Y_2)\in\mathcal{L}(w_2),X_3\in\mathbb{R}\}.$  The plane $\mathcal{PL}(w_2)$ intersects the ellipsoidal manifold $\mathcal{M}(\Escr)$ at a nearly-elliptic closed curve $\mathcal{C}_{PLM}(\Escr,w_2)$, while it intersects the sphere $\mathcal{I}_0$ at the circle $\mathcal{C}_{PLI_0}: s^2+X_3^2=2\AMD$. The curve $\mathcal{C}_{PLM}(\Escr,w_2)$ has central symmetry, and, owing to the fact that its quadratic-form approximation is negative definite, its overall size decreases as $\Escr$ increases from the most negative possible value $\Escr=\Escr_{min}$ towards the value $\Escr=0$, at which $\mathcal{C}_{PLM}(\Escr,w_2)$ reduces to a point at the origin of the plane $\mathcal{PL}(w_2)$. Finally, the surface $\dot{Y}_3=0$ also intersects the plane $\mathcal{PL}(w_2)$ at a curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$. As a result, fixing the value of $w_2$, there are three possibilities as regards the intersections of the curve $\mathcal{C}_{PLM}(\Escr,w_2)$, which varies with the energy, and the curves $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$ and $\mathcal{C}_{PLI_0}$:

Regime 1: the curve $\mathcal{C}_{PLM}(\Escr,w_2)$ intersects the circle $\mathcal{C}_{PLI_0}$ at four points, two of which ($P_1,P_2$) are above the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$, hence corresponding to $\dot{Y}_3>0$, while it intersects the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$ itself at two points which are exterior to the circle $\mathcal{C}_{PLI_0}$ (Fig.<ref>). In this case, the permissible initial conditions in the Poincaré surface of section are given by $X_2=s_2\cos(w_2)$, $Y_2=s_2\sin(w_2)$ with
$s_2$ in the interval $s_{P_1}\leq s_2\leq s_{P_2}$, where $s_{P_1}$, $s_{P_2}$ are the values of $s_2$ to which project the points of intersection $P_1,P_2$. At both these points we have $\dot{Y}_3>0$.

Left: intersections between the surfaces $\mathcal{M}(\Escr)$ at the energy  $\Escr=-9.16\cdot 10^{-5}$ (light gray), the sphere $\mathcal{I}_0$ (thin gray), and the surface $\dot{Y}_3=0$ (green slightly inclined). The thick black curve projected on the plane $(X_2,Y_2)$ defines the permissible domain of initial conditions (see text). Right: the intersections between the surface $\mathcal{M}(\Escr)$ (thick light gray closed curve), the sphere $\mathcal{I}_0$ (thin gray curve and the condition $\dot{Y}_3=0$) (green slightly inclined curve) as viewed in the plane $(X_2,X_3)$ for $Y_2=0$.  The points of intersection mark the limits of possible initial conditions along $X_2$ for $X_3=Y_2=0$.

Regime 2: for larger energies, the curve $\mathcal{C}_{PLM}(\Escr,w_2)$ still intersects the circle $\mathcal{C}_{PLI_0}$ at four points, two of which ($P_1,P_2$) are above the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$, while now intersecting the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$ itself at two points $(P_3,P_4)$, which are interior to the circle $\mathcal{C}_{PLI_0}$ (Fig.<ref>). In this case, the permissible initial conditions in the Poincaré surface of section are given by $X_2=s_2\cos(w_2)$, $Y_2=s_2\sin(w_2)$ with
$s_2$ in one of the intervals $s_{P_3}\leq s_2\leq s_{P_1}$, or $s_{P_2}\leq s_2\leq s_{P_4}$, where $s_{P_1}$, $s_{P_2}$,  $s_{P_3}$, $s_{P_4}$ are the values of $s_2$ to which project the points of intersection $P_1$,$P_2$,$P_3$,$P_4$. In this case we have $\dot{Y}_3>0$ at the inner limits $P_1$, $P_2$, while we have $\dot{Y}_3=0$ (i.e. the orbit arrives tangently to the surface of section) at the outer limits $P_3$, $P_4$.

Same as in Fig.<ref> but for the surface $\mathcal{M}(\Escr)$ calculated at $\Escr=-7.09\cdot 10^{-5}$. The limits of possible initial conditions along $X_2$ for $X_3=Y_2=0$ are marked by the intesections $\mathcal{M}(\Escr)$ both with $\mathcal{I}_0$ and the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$, for $w_2=0$ (see text).

Regime 3: for still larger energies, the curve $\mathcal{C}_{PLM}(\Escr,w_2)$ decreases in size in such a way that it no longer intersects the circle $\mathcal{C}_{PLI_0}$. Then the only limits are posed by its intersections with the curve $\mathcal{C}_{PL\dot{Y}_3=0}(w_2)$ at the points $(P_3,P_4)$, which remain interior to the circle $\mathcal{C}_{PLI_0}$ (Fig.<ref>). Then, the permissible initial conditions in the Poincaré surface of section are given by $X_2=s_2\cos(w_2)$, $Y_2=s_2\sin(w_2)$ with
$s_2$ in the interval $s_{P_3}\leq s_2\leq s_{P_4}$, and both limits correspond to $\dot{Y}_3=0$, i.e. to orbits  tangent to the surface of section.

Inner tangency between the surface $\mathcal{M}(\Escr_{2,3})$ at the energy  $\Escr_{2,3}=-6.77\times 10^{-5}$ and the sphere $\mathcal{I}_0$. The point of tangency marks a fixed point which corresponds to a planar periodic orbit of the type `mode B' (aligned apsidal corotation).

Having fixed the value of the energy $\Escr$, and repeating the computation of the above intersection points for various lines $\mathcal{L}(w_2)$ (i.e. various choices of the argument of the perihelion $\omega_2$), we can explicitly compute the limits of the whole domain $\mathcal{D}(\Escr)$ and proceed in the computation of the phase portraits, obtaining several initial conditions within the domain $\mathcal{D}(\Escr)$. In practice, owing to the symmetries of the Hamiltonian, it is sufficient to consider only a scanning of initial conditions along the line $\mathcal{L}(w_2=0)$, i.e., $X_2=s_2$, $Y_2=0$. Note also that the regimes 1 and 2, as well as 2 and 3, are separated at the energies $\Escr_{1,2}$ and $\Escr_{2,3}$ respectively. We will see in the next section the topological differences in the phase portraits between the various regimes. We only note here that the generic regime is regime 3, which emerges beyond the energy $\Escr_{2,3}$, at which the manifold of constant energy $\mathcal{M}(\Escr_{2,3})$ has an inner tangency with the sphere $\mathcal{I}_0$. Finally, we note that both critical tangencies occuring at the energies $\Escr_{min}$ and $\Escr_{2,3}$ correspond to basic periodic orbits of the system (the apsidal corotation resonances), while the final limit $\Escr=0$ corresponds to the Kozai-Lidov fixed point of the system, at which $\i_{mut}=\i_{max}$ and $\e_2=\e_3=0$.

§.§.§ Poincaré surface of section: precision tests

Figure <ref> shows a summary of our basic example of computation of phase portraits, in the form of Poincaré surfaces of section computed as explained in the previous subsection, referring to a 3D planetary system with mass, periods and $\AMD$ parameters as those estimated for the $\upsilon-$Andromedae system. We adopt the value $\mathcal{G}=4 \pi^2~AU^3/(yr^2 M_\odot)$ for Newton's gravity constant, as well as the mass parameters
$m_0 = 1.31 M_\odot$, $m_2 = 13.98 M_J$, $m_3 = 10.25 M_J$
($M_\odot=$ mass of the Sun, $M_J=0.0009546 M_\odot=$ mass of Jupiter),
and the semi-major axes $a_2 = 0.829$ AU, $a_3= 2.53$ AU (according to Table 13 of [McArthur et al., 2010]).
For what concerns the constants $L_z$ or $\AMD$, we estimate their values adopting for the two planets' eccentricities, inclinations and arguments of periastron and of the nodes the values proposed in Table 1 of [Deitrick et al., 2015], namely,
$\e_2 = 0.2445$, $\e_3 = 0.316$, $\i_2 = 11.347^\circ$, $\i_3 = 25.609^\circ$, $\omega_2 = 247.629^\circ$, $\omega_3 = 252.991^\circ$, $\Omega_2 = 248.181^\circ$, $\Omega_3 = 11.425^\circ$. As emphasized in the introduction, there are great uncertainties in the observation as regards, in particular, the estimates on the planetary eccentricity and inclination vectors. However, starting with parameters as above, we can have a representative value for the system's $\AMD$, which then serves to analyze the secular orbital dynamics under different assumptions for the planets' initial conditions. The value of the $\AMD$ is estimated from the above data, by first computing the angular momentum vector with respect to a heliocentric frame of reference whose x-axis points to the observer, then computing the rotation matrix connecting the observer's frame with Laplace's reference frame, and finally rotating the positions and velocities of both planets to the Laplace frame of reference and re-calculating their inclinations and arguments of periastron and of the nodes. This yields the values $\i_2 = 18.4748^\circ$, $\i_3 = 14.6462^\circ$, $\omega_2 = 289.049^\circ$, $\omega_3 = 235.464^\circ$, $\Omega_2 = 217.318^\circ$, $\Omega_3 = 37.3176^\circ$. Then, we obtain $Lz=0.183101~AU^2 M_\odot/yr$, which, together with $\Lambda_2=0.0873819$, $\Lambda_3=0.111923$, leads to $\AMD=L_2+L_3-L_z=0.0162044~AU^2 M_\odot/yr$.

Poincaré surfaces of section in the plane $(\e_2 \cos(\omega_2 ), \e_2 \sin(\omega_2 ))$ with $L_z$ fixed and different values of energy. The surfaces of section are computed by a numerical integration of trajectories in the Hamiltonian truncated at: Top multipolar degree $N_{P}=5$, order $N_{bk}=10$ in the eccentricities, and energies (from left to right) $\Escr= -6.67\cdot 10^{-5},-2.53\cdot 10^{-5},-1.9\cdot 10^{-5},-1.17\cdot 10^{-5},-2.61\cdot 10^{-6}$. (Middle) $N_{P}=6$, $N_{bk}=10$, and energies top $\Escr= -6.75\cdot 10^{-5}, -2.69\cdot 10^{-5}, -1.9\cdot 10^{-5},-1.16\cdot 10^{-5}, -2.61\cdot 10^{-6}$. (Bottom) $N_{P}=6$, $N_{bk}=12$, and energies $\Escr= -6.75\cdot 10^{-5},-2.69\cdot 10^{-5},-1.9\cdot 10^{-5},-1.16\cdot 10^{-5},-2.59\cdot 10^{-6}$.

In Fig.<ref> we already observe some of the significant changes in the phase portraits when the energy is varied from a value close to $\Escr_{min}\simeq-1.2\cdot 10^{-4}$ (slightly different in each of the models considered in the figure), to another close to $\Escr_{max}=0$. In the next section, we will return to a detailed discussion of the changes observed in the phase portraits as the energy moves from a value, in which the system is closer to the planar regime, to another, in which the system is highly inclined. Here, instead, we only discuss how these figures control the robustness of the phase portraits with respect to the model chosen, which can differ in the maximum multipole degree $N_{P}$ up to which the original Hamiltonian $\Hscr$ is expanded, before the averaging, as well as in the maximum order in book-keeping $N_{bk}$ of the Hamiltonian $\Hscr_{sec}$, which is also equal to the maximum order at which the eccentricities appear in $\Hscr_{sec}$. Checking with increasing values of $N_{P}$ and $N_{bk}$, we find that the phase portraits stabilize at $N_{P}=5$, $N_{bk}=10$. Beyond these values, the sequence of bifurcations of new fixed points and the corresponding changes in the topological features of the phase portraits become marginal, with only changes in the second significant figure observed as regards both the value of the energy $\Escr$ where a bifurcation occurs, and the position of the corresponding fixed points, separatrices, etc.

Since we are interested only in a qualitative description of the structure and evolution of the phase portraits, for reducing computational time we choose as our basic model the one with $N_{P}=5$, $N_{bk}=10$. As an independent test, we perform a comparison between the phase portraits obtained with this model and those obtained by a completely independent Laplace-Lagrange expansion of the Hamiltonian $\Hscr$ averaged over the fast angles. In the latter case, the book-keeping process described in subsection <ref> must be altered for the decomposition $\Hscr=\Hscr_{planar}+\Hscr_{space}$ to naturally emerge while the Jacobi reduction is performed. Namely, we first introduce the book-keeping $\e_j\rightarrow\varepsilon \e_j$, $j=2,3$ and expand the direct term $-\mathcal{G}m_2m_3/|\mathbf{r}_2-\mathbf{r}_3|$ in powers of the orbital eccentricities up to the book-keeping order $N_{bk}=10$. The so-computed expression has the form:

\begin{align}\label{term.direct.part0}
\sum_{\alpha,\beta,\gamma}
\frac{D_{\alpha,\beta,\gamma}(a_2,a_3,\e_2,\e_3,\i_2,\i_3;\varepsilon)}{(\Delta_{2,3})^{\frac{2\,s+1}{2}}} \cos(\alpha_2\lambda_2 + \alpha_3\lambda_3 + \beta_2\Omega_2 + \beta_3\Omega_3 + \gamma_2\varpi_2 + \gamma_3\varpi_3)\, ,
\end{align}

where the denominator $\Delta_{2,3}$, after the substitution $\Omega_3=\Omega_2+\pi\,$, takes the form:

\begin{equation}
\label{den.prima.riduz}
\begin{split}
\Delta_{2,3}=&a_2^2 + a_3^2 - a_2 a_3 \cos(\lambda_2 - \lambda_3) - a_2 a_3  \cos(\lambda_2 - \lambda_3)\cos(\i_2+\i_3)\\
+&a_2 a_3  \cos(\lambda_2 + \lambda_3 - 2 \Omega_2)\cos(\i_2+\i_3) -
a_2 a_3 \cos(\lambda_2 + \lambda_3 - 2 \Omega_2)\, ,
\end{split}\,
\end{equation}

with $\varpi_j=\omega_j\,+\,\Omega_j$ and $\lambda_j=M_j\,+\,\varpi_j$ denoting the longitudes of the pericenter and the mean longitudes of the bodies $j=2,3$ respectively. The integer vectors $\alpha\equiv(\alpha_2,\alpha_3)$, $\beta\equiv(\beta_2,\beta_3)$, $\gamma\equiv(\gamma_2,\gamma_3)$ satisfy the D'Alembert rule $\alpha_{2}+\alpha_{3}+\beta_2+\beta_3+\gamma_2+\gamma_3=0\,$. The coefficients $D_{\alpha,\beta,\gamma}$ are polynomial in the book-keeping parameter $\varepsilon$. At this point, we introduce the book-keeping identities:

\cos(i_2+i_3)=\eta\varepsilon^2(\cos(i_2+i_3)-1)+1\, ,

i.e. (see Eq (<ref>))

\begin{align*}
&\cos(\i_2)\cos(\i_3)=\varepsilon^2\eta\left(\cos(\i_2)\cos(\i_3)-1\right)+1\, , &
&\sin(\i_2)\sin(\i_3)=\varepsilon^2\eta\sin(\i_2)\sin(\i_3)\, ,
\end{align*}

and expand again the expression (<ref>) with respect to the book-keeping parameter $\varepsilon$ up to the truncation order $N_{bk}=10$. This leads to an expression of the form

\begin{align}\label{term.direct.part}
\frac{C(a_2,a_3,\e_2,\e_3,\i_2,\i_3;\varepsilon,\eta)}{\left(a_2^2 + a_3^2 - 2 \,a_2\,a_3\, \cos(\lambda_2 - \lambda_3)\right)^{\frac{2\,s+1}{2}}} \cos(\alpha_2\lambda_2 + \alpha_3\lambda_3 + (\beta_2+\beta_3)\Omega_2  + \gamma_2\varpi_2 + \gamma_3\varpi_3)\, ,
\end{align}

where the coefficients $C(a_2,a_3,\e_2,\e_3,\i_2,\i_3;\varepsilon,\eta)$ are polynomial in the book-keeping parameters $\varepsilon,\eta$. Finally, we perform the classical Laplace-Lagrange averaging `by scissors':
the denominator of (<ref>) is Fourier-expanded

\begin{equation}
\label{espans.flam}
\frac{1}{\left(a_2^2 + a_3^2 - 2 \,a_2\,a_3\, \cos(\lambda_2 - \lambda_3)\right)^{\frac{2\,s+1}{2}}}=a_3^{-(2\,s+1)}\sum _{j\geq 0} b_{s+\frac{1}{2}}^{(j)}\left(\frac{a_2}{a_3}\right)\cos(j(\lambda_2-\lambda_3))\, ,
\end{equation}

where $b_{s+\frac{1}{2}}^{(j)}$ are Laplace coefficients, which can be computed numerically, as

\begin{equation}
\begin{split}
\label{Lapl.coeff}
&b_{s+\frac{1}{2}}^{(0)}\left(\alpha\right)=\frac{1}{2\pi}\int_{0}^{2\pi} \left(1+\alpha^2 - 2 \,\alpha\, \cos(\vartheta)\right)^{-(s+\frac{1}{2})}\,d\vartheta\, ,\\
&b_{s+\frac{1}{2}}^{(j)}\left(\alpha\right)=\frac{1}{\pi}\int_{0}^{2\pi} \left(1+\alpha^2 - 2 \,\alpha\, \cos(\vartheta)\right)^{-(s+\frac{1}{2})}\cos(j\vartheta)\,d\vartheta \qquad j\geq 1\, ,
\end{split}
\end{equation}

where $\alpha=a_2/a_3\,$, or via a multipolar expansion (see the Appendix <ref>).
After performing the above expansions, the secular Hamiltonian is obtained by dropping from the original Hamiltonian all terms depending on the `fast' angles $\lambda_i$. This leads to a model in which the inclinations appear only through the combination $\cos(i_2+i_3)$. Switching, now, back to the angles $\omega_j,\Omega_j$, performing the Jacobi reduction as in subsection <ref> (i.e., using Eq. (<ref>), collecting together all terms independent or depending on $\eta$, and restoring the numerical values of the book-keeping coefficients $\varepsilon=\eta=1\,$), we arrive at a `Laplace-Lagrange' model $\Hscr_{LL}=\Hscr_{LL,planar}+\Hscr_{LL,space}$ which has the same form as the model of Eq. (<ref>). This is further processed with the introduction of Poincaré variables as in step 4 of subsection <ref>.

Comparison of the Poincaré surfaces of section in the plane $(\e_2 \cos(\omega_2 ), \e_2 \sin(\omega_2 ))$ with $L_z$ fixed and different values of the energy, in two models: Top: our basic model $N_{P}=5$, $N_{bk}=10$ (same as in the top row of Fig.<ref>). Bottom: the Laplace-Lagrange secular Hamiltonian model $\Hscr_{LL}$ (see text) truncated at order 10 in the eccentricities, with energies (from left to right) $\Escr= -6.62\cdot 10^{-5}, -2.94\cdot 10^{-5}, -1.92\cdot 10^{-5},-1.18\cdot 10^{-5}, -2.73\cdot 10^{-6}$.

From Fig.<ref> we conclude that similar remarks as those of Fig.<ref> can be made as regards the comparison of the secular model $\Hscr_{sec}$ adopted in the present work, based on a multipolar expansion, and the model $\Hscr_{LL}$ obtained by the classical Laplace-Lagrange expansion in the eccentricities. We note, however, that the latter is much harder to compute, while it is obtained by using a series reversion of Kepler's equation which has a limited convergence. On the other hand, the model based on closed-form averaging of the multipole Hamiltonian expansion requires a special treatment as regards the computation of the underlying canonical transformation which eliminates (instead of `scissor cutting') the fast angles (see [Cavallari and Efthymiopoulos, 2022]). Since in the present work we do not consider this transformation, we will hereafter deal with results derived only by use of the simple multipole and closed-form averaged model $\Hscr_{sec}$, with the truncation orders $N_{P}=5$, $N_{bk}=10$.

§ DYNAMICS

§.§ General

In this section we are interested in analyzing the most important phenomena observed in the phase portraits of the Hamiltonian $\Hscr_{sec}(X_2, X_3, Y_2, Y_3; \AMD)$ computed with the basic reference model corresponding to the truncation orders $N_{P}=5$, $N_{bk}=10$. In the following subsections we will first present the general picture of the transitions taking place in the structure of the phase space as the energy increases in the range $\Escr_{min}\leq\Escr\leq 0$. Such transitions are caused, for example, by changes in the nature (e.g. stability) of the main equilibria of the system, giving birth to new families of periodic orbits which connect the families dominant in the planar-like regime (high eccentricities, low mutual inclination) with those of the highly inclined regime (low eccentricities, high mutual inclination).

In the analysis of phase portraits as in the sequel, we parametrize all transitions due to bifurcations of new periodic orbits using the fixed (in the Poincaré section) energy $\Escr$ as the parameter. All energies referred to below are given in units of $M_\odot AU^2/yr^2$. It is easy to see that a certain value of the energy $\Escr$ establishes a range of allowed mutual inclinations

\begin{equation}
\i_{mut}^{min}(\Escr)\leq \i_{mut}\leq \i_{mut}^{max}(\Escr)
\end{equation}

where both $\i_{mut}^{min}(\Escr)$ and $\i_{mut}^{max}(\Escr)$ are increasing functions of $\Escr$, as shown in Fig.<ref>.

The minimum and maximum possible values of the mutual inclination $\i_{mut}^{min}$, $\i_{mut}^{max}$, as a function of the energy $\Escr$. The vertical lines correspond to the energies $\Escr_{min}=-1.18237\cdot 10^{-4}$ (birth of the mode A - anti-aligned apsidal corotation), $\Escr_{1,2}=-8\cdot 10^{-5}$ (completion of the libration domain around the mode A), $\Escr_{2,3}=-6.77\cdot 10^{-5}$ (birth of the mode B - aligned apsidal corotation), $\Escr_C=-2.7\times 10^{-5}$ (bifurcation of the inclined Kozai-Lidov periodic orbits $C_1,C_2\,$), and $\Escr_{C,2}=-5.0\cdot 10^{-6}$ (inclined orbit $C_2$ becomes unstable, see text).

Phase portraits (Poincaré surfaces of section in the plane $(\e_2 \cos(\omega_2 ), \e_2 \sin(\omega_2 ))$) in the basic model $N_{P}=5$, $N_{bk}=10$, and for the energies (from from top left to bottom right) $\Escr= -1.145\cdot 10^{-4}$, $-9.16\cdot 10^{-5}$, $-8.23\cdot 10^{-5}$, $-7.09\cdot 10^{-5}$,
$-6.77\cdot 10^{-5}, -5.72\cdot 10^{-5}, -4.93\cdot 10^{-5},-3.81\cdot 10^{-5},-3.19\cdot 10^{-5},-2.76\cdot 10^{-5},-2.53\cdot 10^{-5},-2.16\cdot 10^{-5}, -2.08\cdot 10^{-5}, -1.9\cdot 10^{-5}, -1.58\cdot 10^{-5}, -1.53\cdot 10^{-5}, -1.17\cdot 10^{-5}, -7.69\cdot 10^{-6}, -2.61\cdot 10^{-6}, -7.39\cdot 10^{-7}$. The positions of the fixed points corresponding to the periodic orbits of the modes $A$, $B$, $C_1$, $C_2$, $D_1$, $D_2$ (see text) are marked by arrows.

The minimum and maximum possible values of the mutual inclination $\i_{mut}^{min}$, $\i_{mut}^{max}$, as a function of the energy $\Escr$ are computed as follows: consider the ellipsoidal surface $\mathcal{I}_{\i_{mut}}$ defined by setting $Y_3=0$ in Eq.(<ref>) for a certain value of $\i_{mut}$. Consider the critical energy $\Escr=\Escr_{2,3}$, equal to $\Escr_{2,3}=-6.77\times 10^{-5}$ in our example. For energies $\Escr\geq\Escr_{2,3}$, we compute the values $\i_{mut}=\i_{mut}^{min}(\Escr)$ and $\i_{mut}=\i_{mut}^{max}(\Escr)$ for which the ellipsoid $\mathcal{I}_{\i_{mut}}$ comes tangent to the ellipsoidal manifold of constant energy $\mathcal{M}(\Escr)$, with $\mathcal{I}_{\i_{mut}}$ being at the interior of $\mathcal{M}(\Escr)$ (for $\i_{mut}=\i_{mut}^{max}$) or at the exterior of $\mathcal{M}(\Escr)$ (for $\i_{mut}=\i_{mut}^{min}$). For energies $\Escr_{min}\leq\Escr<\Escr_{2,3}$ (with $\Escr_{min}=-1.182\times 10^{-4}$ in our example), the condition of tangency of $\mathcal{I}_{\i_{mut}}$ at the exterior of $\mathcal{M}(\Escr)$ leads to unphysical values $\cos(\i_{mut})>1$. Thus, in this interval of energies we have $\i_{mut}^{min}=0$, while $\i_{mut}^{max}$ keeps being represented by an increasing function of $\Escr$ (Fig.<ref>).

The vertical lines in Fig. <ref> indicate values of the energy where important changes take place in the structure of the phase portraits due to the birth, or change of stability character, of some of the most important families of periodic orbits of the system. The most important transitions taking place in the structure of the phase portraits are shown in Fig.<ref>, whose details will be presented in subsequent subsections. Increasing the energy $\Escr$, these transitions appear in summary by the following sequence:

i) At the energy $\Escr=\Escr_{min}$, the available domain $\mathcal{D}(\Escr)$ reduces to a point, corresponding to the point of tangency of Fig.<ref>. This is a fixed point of the Poincaré map, whose associated orbit yields two coplanar ellipses with anti-aligned pericenters precessing by the same frequency, known as the apsidal corotation orbit (see [Laughlin et al., 2002], [Lee and Peale, 2003], [Beaugé et al., 2003]). For energies $\Escr>\Escr_{min}$, the above fixed point is continued by a family of periodic orbits, called below the mode A. These correspond physically to inclined planetary orbits whose eccentricities undergo small periodic oscillations around some non-zero constant values $\e_{2,A}$, $\e_{3,A}$ (functions of the energy), while the arguments of perihelia undergo small periodic oscillations around the fixed relation $\omega_2-\omega_3=0$ (see subsections <ref> and <ref>). For energies $\Escr_{min}<\Escr\leq\Escr_{2,3}$ (equal to $-6.77\cdot 10^{-5}$ in our numerical example), the mode A, which generalizes the anti-aligned apsidal corotation family to the non-planar case, is the unique important stable family in the surface of section (see first three panels of top row of Fig.<ref>). The corresponding fixed point is surrounded by closed invariant curves, which represent orbits performing quasi-periodic oscillations around the configuration of anti-aligned perihelia. Up to the energy $\Escr_{1,2}$ (equal to about $-8\cdot 10^{-5}$ in our example), which marks the transition from two to four limits of permissible motion as in Figures <ref> and <ref>, only quasi-periodic orbits around the A mode exist. On the other hand for energies $\Escr_{1,2}<\Escr<\Escr_{2,3}$, we can also have trajectories with argument $\omega_3-\omega_2$ either circulating or librating around the value $\omega_3-\omega_2=\pi$ (alignment). As explained in detail in subsection <ref>, the separation between the various librating or circulating regimes is not due to the presence of a dynamical separatrix, but can be explained by an integrable Hamiltonian model approximating $\Hscr_{sec}$ in the corresponding energy regime, whose phase space has the topology of a 3-sphere rather than the plane $\mathbb{R}^2$.

ii) For energies $\Escr_{1,2}<\Escr<\Escr_{2,3}$ there is a prohibited domain surrounding the center of the librating motions around $\omega_3-\omega_2=\pi$ (see fourth panel, top row of Fig.<ref>). At the energy $\Escr=\Escr_{2,3}$ this domain shrinks to zero, and at the center of the librations appears a second fixed point of the Poincaré map, which is a periodic orbit physically corresponding to the planar aligned apsidal corotation orbit (see the first panel in the second row of Fig. <ref>). This also marks the inner point of tangency of the sphere $\mathcal{I}_0$ with the energy manifold $\mathcal{M}(\Escr_{2,3})$ (Fig.<ref>). As indicated in Fig.<ref>, for energies $\Escr>\Escr_{2,3}$ there can no longer be any planar orbit intersecting the surface of section. However, similarly as for the mode A, the fixed point corresponding to the aligned apsidal corotation is continued as a family of off-plane periodic orbits, hereafter called the mode B. Physically, such orbits undergo small periodic oscillations around some non-zero constant values $\e_{2,B}$, $\e_{3,B}$ (also being both functions of the energy), while the arguments of perihelia undergo small periodic oscillations around the fixed relation $\omega_2-\omega_3=\pi$ (see also subsections <ref>, <ref>). The mode B also is surrounded by quasi-periodic orbits with arguments of the perihelia librating around the relation $\omega_3-\omega_2=\pi$.

iii) The topology of the phase space induced by the alteration between the domains of libration and circulation around the modes A and B dominates the picture obtained for the phase portraits in a large subinterval within the permissible range of values of the energy (up to about $\Escr=-3\cdot 10^{-5}$ in our example). As seen in Fig.<ref>, this covers most cases of highly-inclined orbits, with mutual inclinations being in our example as high as $\sim 40^\circ$. However, at a critical energy $\Escr_{C}$ (equal to about $-2.7\cdot 10^{-5}   $ in our example), a saddle-node bifurcation takes place, giving rise to two new fixed points of the Poincaré map, corresponding to periodic orbits hereafter called the Kozai-Lidov orbits $C_1$ and $C_2$. Physically, these are highly inclined orbits with planetary eccentricities $\e_2$ and $\e_3$ smaller than those of the modes A and B, tending actually to zero as $\Escr\rightarrow 0$. The detailed sequence of bifurcations related to these orbits is discussed in detail in subsection <ref> (see panels 11 to 20 of Fig.<ref>). The most important transitions regard the orbit $C_1$, which becomes stable nearly immediately after its birth, while the orbit $C_2$ undergoes the classical Lidov-Kozai transition from stability to instability, accompanied by the appearance of chaotic motions around it. Such phenomena appear in a small range of energies near the limit $\Escr=0$, where the whole phase space shrinks again to a unique point at the origin of the surface of section, corresponding to two circular orbits having the maximum possible mutual inclination ($\i_{mut}^{max}(\Escr=0)\simeq 46^\circ\,$ in our example). In fact, Fig.<ref> indicates that in this Lidov-Kozai regime, the range in possible inclinations for all orbits becomes narrow, being limited to values around the critical $\Escr_{C,2}$ where the periodic orbit $C_2$ underhgoes the Lidov-Kozai instability.

The above was a synoptic account of the main phenomena occuring in the phase space as the energy (and, hence, mutual inclination) of the orbits increases and the system gradually shifts from a regime where the dynamics is qualitatively similar to the one of the planar problem, to the Lidov-Kozai regime dominated by nearly circular and highly inclined orbits. In the next subsections we examine these phenomena in more detail, following the same sequence as the one of their appearance by increasing the value of the energy $\Escr$.

§.§ Planar-like regime

In the present section we will discuss in detail the structure of the phase portraits in the range of energies $\Escr_{2,3}\leq\Escr<\Escr_C$, were the dominant periodic orbits are the Modes A and B, which generalize the apsidal corotations (anti-aligned and aligned, respectively) of the planar case. It was already mentioned that the phase portait in this case contains two domains where the argument $\omega_2-\omega_3$ librates (around the values $0$ and $\pi$ respectively), separated by a domain where $\omega_2-\omega_3$ circulates, as in panels 5 to 10 of Fig.<ref>. Owing to its similarity with the phase portrait of the planar problem, this will be called the planar-like regime. It was emphasized, however, that the maximum mutual inclination can be quite high in this regime (see Fig.<ref>), thus the analogy with the planar case stems from the dynamics, and not necessarily from the degree of coplanarity of the orbits in this regime. Excluding the energetically prohibited domains (as specified in subsection <ref>), this can be extended to cover the cases where the librational domains are only partly covered by quasi-periodic orbits, as in panels 1-4 of Fig.<ref>.

§.§.§ Integrable approximation of the Hamiltonian

The main qualitative features of the planar-like regime, as well as precise computations regarding its periodic and surrounding quasi-periodic orbits, can be obtained in the context of an integrable approximation for the Hamiltonian $\Hscr_{sec}$, stemming from the splitting $\Hscr_{sec}=\Hscr_{planar}+\Hscr_{space}$ as in Eq. (<ref>). Starting from

\begin{equation}
\Hscr_{sec}(a_2, a_3, \e_2, \e_3, w_2, w_3; L_z)\,=\Hscr_{planar}(w_2-w_3,W_2,W_3)
+\Hscr_{space}(w_2,w_3,W_2,W_3; L_z)\, .
\end{equation}

and splitting $\Hscr_{space}$ into those terms which depend only on the difference $w_2-w_3=\omega_3-\omega_2$, denoted by $\Hscr_{0,{space}}$, and those which do not, denoted by $\Hscr_{1,{space}}$, we arrive at the following decomposition of the Hamiltonian

\begin{equation}\label{hamdecompo}
\Hscr_{sec}=\underbrace{\Hscr_{planar}( w_2-w_3,W_2,W_3) + \Hscr_{0,space}(w_2-w_3,W_2,W_3; L_z)}_{\textbf{\textit{ Integrable part}}:=\Hscr_{int}}+ \Hscr_{1,space}(w_2,w_3,W_2,W_3; L_z)\, .
\end{equation}

The first two terms in the above expression give rise to a 2 degrees of freedom integrable Hamiltonian

\begin{equation}\label{hamint}
\Hscr_{int}( w_2-w_3,W_2,W_3;\AMD)= \Hscr_{planar}( w_2-w_3,W_2,W_3)+ \Hscr_{0,space}( w_2-w_3,W_2,W_3;\AMD)
\end{equation}

whose second integral is $W_2+W_3$.

A fact hidden in the process of Jacobi reduction is that a decomposition of the Hamiltonian as in Eq.(<ref>) yields a relative importance of the terms $\Hscr_{int}$ and $\Hscr_{1,space}$ varying
with the energy level $\Escr$ at which the orbits are computed. This is due to the fact that, after throwing apart constants, all the terms in $\Hscr_{space}$ stem from substitutions of the inclinations $\i_2$ and $\i_3$ depending only on the small quantity $1-\cos(\i_2+\i_3)$ and being of order second or higher in the eccentricities, according to Eq. (<ref>) and (<ref>). Thus, both $\Hscr_{0,space}$ and $\Hscr_{1,space}$ contain terms weighted by factors $(1-\cos(\i_2+\i_3))^s(\e_j^2+\ldots)$ (with $s\geq 1$, $j=2,3$). On the other hand, the terms $\Hscr_{planar}$ contain no factors $(1-\cos(\i_2+\i_3))^s$ and are of degree quadratic or higher in the eccentricities. Thus, due to the bound between increasing energy $\Escr$ and increasing mutual inclination (Fig.<ref>), the relative importance of the terms $\Hscr_{1,space}$ with respect to the terms $\Hscr_{int}$ in the Hamiltonian rises as the energy (and hence the level of mutual inclination) increases. This is demonstrated graphically in Fig.<ref>, which shows a comparison between the form of the manifolds
of constant energy $\mathcal{M}(\Escr)$ and $\mathcal{M}_{int}(\Escr)$ (Eqs.(<ref>) and (<ref>)) computed at four different energy levels chosen as $\Escr_{min}<\Escr_1<\Escr_{2,3}$, $\Escr_{2,3}<\Escr_2<\Escr_{C,2}$, $\Escr_3=\Escr_{C,2}$, $\Escr_{C,2}<\Escr_4<0$. While the manifolds of $\Hscr_{int}$ remain always ellipsoidal-like, we note the progressive change of the form of the manifolds of constant energy in the complete model from an ellipsoidal to a peanut-shaped form, as the energy increases. This is caused by the growing importance of some terms in $\Hscr_{1,space}$, in particular the terms $\cos(2\omega_2)$ quadratic in $\e_2$, which are the same terms causing the transition to the Lidov-Kozai regime (see subsection <ref>). In fact, near the energy $\Escr_{C,2}$, where the orbit $C_2$ turns from stable to unstable, the manifold $\mathcal{M}(\Escr)$ under the complete model becomes nearly cylindrical, marking the change of its section with the plane $(X_2,Y_2)$ from elliptic-like to hyperbolic-like, as implied by the Lidov-Kozai mechanism.

The top row shows the manifolds of constant energy $\mathcal{M}(\Escr)$ compared to the manifolds of constant energy $\mathcal{M}_{int}(\Escr)$ of the integrable model $\Hscr_{int}$ (bottom row), for the energies  (from left to right) $\Escr_1=-8\cdot 10^{-5}$ ($\Escr_1<\Escr_{2,3}$), $\Escr_2=-4\cdot 10^{-5}$ ($\Escr_{2,3}<\Escr_2<\Escr_{C,2}$), $\Escr_3=-5\cdot 10^{-6}$ ($\Escr_3=\Escr_{C,2}$), $\Escr_4=-5\cdot 10^{-7}$ ($\Escr_4>\Escr_{C,2}$).

Given the above, we will now focus on a description of the phase portraits in the energy regime (roughly identified as $\Escr<\Escr_C$) where the dynamics induced by $\Hscr_{sec}$ can be well approximated by the dynamics of $\Hscr_{int}$. In this regime, the following canonical transformation proves useful in semi-analytical (normal form) calculations related to the periodic orbits A and B and their surrounding quasi-periodic orbits:

\begin{align}
\label{trasf.can.int.plane}
& \psi=w_2-w_3\, , & &\Gamma=\frac{W_2-W_3}{2}\, , \notag\\
& \varphi=w_2+w_3\, , & &J=\frac{W_2+W_3}{2}~~.
\end{align}

The Hamiltonian in the new variables reads (apart from a constant)

\begin{equation}
\label{Ham.new.var}
\Hscr_{sec}(\psi, \varphi,\Gamma, J)=\Hscr_{int}( \psi, \Gamma; J)+\eta\,\Hscr_{1,space}(   \psi, \varphi, \Gamma, J)\, .
\end{equation}

Figure <ref> shows the phase portrait (in the representative planes $(\e_2 \cos(\omega_2), \e_2 \sin(\omega_2)$ and $(\omega_2, \e_2)\,$) corresponding to the Hamiltonian $\Hscr_{int}$ at the energy $\Escr=-6.6\cdot 10^{-5}$. The phase portrait, computed as a Poincaré surface of section $\mathcal{P}(\Escr=-6.6\cdot 10^{-5}, \, \AMD=0.0162044)$ (see Eq. (<ref>)), yields invariant curves equivalent to those obtained by the continuous flow after treating $\Hscr_{int}$ as a one degree of freedom Hamiltonian in the variables $(\psi,\Gamma)$, with $J$ serving as parameter.

We observe that with the integrable model $\Hscr_{int}$ we obtain a phase portrait with features qualitatively very similar to those of the phase portraits in the `planar-like' regime under the complete Hamiltonian (e.g. the panels 5-10 in Fig.<ref>). In particular, the modes A and B of the integrable model are found as fixed points of $\Hscr_{int}$, given by the solutions of the equations

\begin{equation}
\label{eq.moto.apsidal}
\begin{cases}
\displaystyle{\dot{\psi}=\frac{\partial\Hscr_{int}}{\partial\Gamma}=0}
\\[2ex]
\displaystyle{\dot{\Gamma}=-\frac{\partial\Hscr_{int}}{\partial\psi}=0}
\end{cases} \, .
\end{equation}

Since $J$ is an integral of motion of $\Hscr_{int}$, setting $J=K=$ constant implies that any solution $( \psi_{\ast},\,\Gamma_{\ast})$ satisfying (<ref>) is a periodic orbit with period given by $T_{\varphi}=2\pi/\omega_{\varphi}\,$, where $\displaystyle{\dot{\varphi}=\omega_{\varphi}=\frac{\partial\Hscr_{int}}{\partial J}\vert_{(\psi=\psi_{\ast},\,\Gamma=\Gamma_{\ast}\, , J=K})}\,$. In particular, the modes $A\,$ and $B\,$ are given, respectively, by $(\psi^{(A)}=0,\,\Gamma^{(A)}\, ,K^{(A)})$  and $( \psi^{(B)}=\pi,\,\Gamma^{(B)},\, K^{(B)})\,$, where $(\psi^{(A)},\,\Gamma^{(A)})\,$, $(\psi^{(B)},\,\Gamma^{(B)})\,$ are solutions of (<ref>). In fact, since the Hamiltonian $\Hscr_{int}$ depends only on the harmonics $\cos(k\psi)$, with $k\in\naturali\,$, the partial derivative in the second of Eqs. (<ref>) yields only $\sin(k\psi)\,$ terms, hence, it admits the solution $\psi=\psi^{(A)}=0$ and $\psi=\psi^{(B)}=\pi$. Given now that $\Omega_3=\Omega_2+\pi$, the condition $\omega_3=\omega_2\,$, i.e. $\psi=0=\psi^{(A)}\,$, implies $\varpi_3=\varpi_2+\pi$ (perihelia anti-aligned), while the condition $\omega_3=\omega_2+\pi\,$, i.e. $\psi=\pi=\psi^{(B)}\,$, implies (modulo $2\pi$) $\varpi_3=\varpi_2$ (perihelia aligned). On the other hand, substituting one of the angles, $\psi^{(A)}$ or $\psi^{(B)}$, in the first of the Eqs. (<ref>), we obtain an algebraic equation of the form $\dot{\psi}(\psi=\psi^{(A)},\Gamma,J)=0$ or $\dot{\psi}(\psi=\psi^{(B)},\Gamma,J)=0$. This, together with the constant energy condition $\Hscr_{int}(\psi=\psi^{(A)},\Gamma,J)=\Escr$, or $\Hscr_{int}(\psi=\psi^{(B)},\Gamma,J)=\Escr$, can be solved to yield the pairs of values $\Gamma^{(A)},J=K^{(A)}$, or $\Gamma^{(B)},J=K^{(B)}$. Finally, the frequency of the apsidal precession for anyone of the periodic orbits A, B is given by $\nu_{apsidal}=-\omega_\varphi/2=-\dot{\varphi}/2=-(1/2)(\partial\Hscr_{int}(\psi,J,\Gamma)/\partial J)$, with $(\psi,J,\Gamma)$ substituted with one of the solutions A or B.

§.§.§ The phase space of $\Hscr_{int}\,$: Hopf variables

An alternative method to compute the equilibria A and B stems from the use of a particular set of variables, called the Hopf variables [Cushman and Bates, 1997], which, besides the computation of the equilibria, provides a global mapping of the phase space of the integrable Hamiltonian $\Hscr_{int}$ to the 3-sphere, thus allowing for a clear identification of all possible orbital dynamical regimes. We introduce the variables $(\sigma_1,\, \sigma_2,\, \sigma_3\,)$ defined by:

\begin{equation}\label{def.Weyl}
\sigma_1=X_2 X_3+Y_2 Y_3\, ,\qquad
\sigma_2=Y_2 X_3-Y_3 X_2\, ,\qquad
\sigma_3=\frac{1}{2}\left(X_2^2+Y_2^2-X_3^2-Y_3^2\right)\, ,
\end{equation}

satisfying the Poisson algebra $\poisson {\sigma_i} {\sigma_j} =-2\,\epsilon_{i j k}\sigma_k\, $, where $\epsilon_{i j k}$ is the Levi-Civita symbol and $i,j,k=1,2,3$. Furthermore, we introduce the variable

\begin{equation}\label{def.sigma0}
\sigma_0=\frac{1}{2}\left(X_2^2+Y_2^2+X_3^2+Y_3^2\right)
\end{equation}

which is a Casimir invariant of the previous algebra, since all Poisson brackets $\{\sigma_i,\sigma_0\}$, $i=1,2,3\,$, vanish. From the definition (<ref>) it follows that

\begin{equation}\label{def.Weyl.Delau}
\sigma_1=2\sqrt{J+\Gamma}\sqrt{J-\Gamma}\cos(\psi)\, ,\qquad
\sigma_2=-2\sqrt{J+\Gamma}\sqrt{J-\Gamma}\sin(\psi)\, ,\qquad
\sigma_3=W_2-W_3=2\Gamma\, .
\end{equation}

We also have the relation $\sigma_0=W_2+W_3=2J$, as well as

\begin{equation}\label{cond.sphere}
\sigma_1^2+\sigma_2^2+\sigma_3^2=\sigma_0^2=4J^2\, .
\end{equation}

Then, given the values of $(\sigma_1,\sigma_2,\sigma_3)$, the values of $\Gamma$, $J$ and $\psi$ can be computed unequivocally using the relations (<ref>) and (<ref>). Furthermore, since $J=\sigma_0/2$, and since the only trigonometric terms in the Hamiltonian $\Hscr_{int}$ are terms $\cos(k(w_3-w_2))=\cos(k\psi)$, $k=\mathbb{N}^{\ast}$, it follows that $\Hscr_{int}=\Hscr_{int}(\sigma_1,\sigma_3;\sigma_0)$, i.e. the Hamiltonian $\Hscr_{int}$ does not depend on $\sigma_2$. This implies that, fixing a value of $\sigma_0$ (i.e. of the integral $J$), the continuous in time phase flow obtained by solving the equations

\begin{equation}\label{sigmaflow}
\dot{\sigma}_1=\{\sigma_1,\sigma_3\}{\partial\Hscr_{int}\over\partial\sigma_3},~~~
\dot{\sigma}_2=
\{\sigma_2,\sigma_1\}{\partial\Hscr_{int}\over\partial\sigma_1}+
\{\sigma_2,\sigma_3\}{\partial\Hscr_{int}\over\partial\sigma_3},~~~
\dot{\sigma}_3=-\{\sigma_1,\sigma_3\}{\partial\Hscr_{int}\over\partial\sigma_1},~~~
\end{equation}

yields a flow equivalent to the one obtained under the Hamiltonian $\Hscr_{int}(\psi,\Gamma;J)$, i.e., treating $J$ as a parameter. Due to the constrain on the sphere (Eq.(<ref>)), the curves of the flow (<ref>) are given by the intesection of the constant energy surface $\Hscr_{int}(\sigma_1,\sigma_3;\sigma_0)=\Escr$ with the sphere, i.e., they are closed curves which can be mapped to invariant curves in the plane $(\Gamma,\psi)$. These are geometrically equivalent to the invariant curves of the Poincaré surface of section of $\Hscr_{int}$ treated as a 2DOF system, mapping $(X_2,Y_2)$ as $X_2=-\sqrt{2(\Gamma+J)}\cos(\psi-\pi)$, $Y_2=\sqrt{2(\Gamma+J)}\sin(\psi-\pi)$. It should be stressed, however, that the above Hamiltonian reduction only yields a geometric equivalence of the two curves, since in the 1DOF reduced system the flow is continuous, while in the 2DOF full system the curves are traced stroboscopically, and they correspond to the intersection of the 2D invariant tori involving both the angles $(\psi,\varphi)$ with the selected surface of section.

Left column: The phase portrait of the integrable Hamiltonian $\Hscr_{int}$ at the energy $\Escr=-6.6\cdot 10^{-5}$, projected in the variables $(\e_2\cos\omega_2,\e_2\sin\omega_2)$ (top), or simply $(\omega_2,\e_2)$ (bottom). Right column, top: At the values $\sigma_0=\sigma_0^{(A)}$ and $\sigma_0=\sigma_0^{(B)}$, the corresponding spheres $\mathcal{S}_{\sigma_0^{(A)}}$, $\mathcal{S}_{\sigma_0^{(B)}}$ become tangent to the energy surfaces $\mathcal{C}_{\sigma_0^{(A)},\,\Escr}$, $\mathcal{C}_{\sigma_0^{(B)},\,\Escr}$. The points of tangency yield the position of the fixed points $A$ and $B$ in the surface of section (see text). Right column, bottom: the intersection of the spheres $\mathcal{S}_{\sigma_0}$ and of the energy surfaces $\mathcal{C}_{\sigma_0,\,\Escr}$ with the plane $(\sigma_1,\sigma_3)$ for $\sigma_2=0$, for various values of $\sigma_0$. The intersection of one sphere with one energy surface yields a curve on the sphere which is projected to a curve in the above plane. For a particular value of $\sigma_0=\sigma_0^{(S)}$, the curve (thick purple) passes through the south pole $S$ of the corresponding sphere $\mathcal{S}_{\sigma_0}$. This corresponds to a trajectory forming a closed curve in the Poincaré section, which surrounds mode A and passes through the origin. This curve delimits the domain of orbits whose angle $\psi=\omega_3-\omega_2$ librates around the value $\psi=0$. At a different value of $\sigma_0=\sigma_0^{(N)}$ the curve of constant energy (thick black in the right column, bottom) passes through the corresponding sphere's north pole $N$. This corresponds to a curve in the surface of section which surrounds the previous curve as well as the origin. In particular, the whole black curve, except for the point $N\,$, yields the part of the corresponding curve in the Poincaré section (left, top figure) contained in the positive semi-plane, while the point $N$ itself inflates to the part of the corresponding closed curve in contained in the negative semi-plane. The domain in the surface of section between the thick purple and the thick black curves corresponds to orbits whose argument $\psi=\omega_3-\omega_2$ circulates. All trajectories beyond the outer delimiting curve exhibit librations of the argument $\omega_3-\omega_2$ around the value $\pi$, characteristic of the B-mode.

The above information can now be used to the purpose of analyzing the structure of the phase portraits using the mapping of each invariant curve in the sphere (<ref>) to the corresponding curve in the usual Poincaré surface of section. To this end, let

\begin{equation}
\label{sphere}
\mathcal{S}_{\sigma_0}=\lbrace(\sigma_1, \sigma_2, \sigma_3)\in\reali^3\, : \, \sigma_1^2+\sigma_2^2+\sigma_3^2=\sigma_0^2 \rbrace
\end{equation}

denote the sphere corresponding to the value $\sigma_0=2J$ of the integral $J$, and

\begin{equation}
\label{energy.curve}
\mathcal{C}_{\sigma_0,\, \Escr} =\lbrace(\sigma_1, \sigma_2, \sigma_3)\in\reali^3\, : \, \Hscr_{int}(\sigma_0, \sigma_1, \sigma_3)=\Escr\rbrace\, ,
\end{equation}

denote the energy surface in the space $(\sigma_1,\sigma_2,\sigma_3)\in\mathbb{R}^3$ corresponding to a fixed energy value $\Escr$. We can have a physical trajectory for all values of $\sigma_0$ (i.e., of $J$) for which the surfaces $\mathcal{S}_{\sigma_0}$ and $\mathcal{C}_{\sigma_0,\, \Escr}$ intersect, limited by two values of $\sigma_0\,$, that are $\sigma_0^{(A)}$ and $\sigma_0^{(B)}$ (corresponding to $J^{(A)}=\sigma_0^{(A)}/2$ and $J^{(B)}=\sigma_0^{(B)}/2\,$), where the two surfaces become tangent (Fig.<ref>). By the non-dependence of $\Hscr_{int}$ on $\sigma_2$, the constant energy surface $\mathcal{C}_{\sigma_0,\, \Escr}$ is normal to any plane $(\sigma_1,\sigma_3)$ with $\sigma_2=const$. Hence, at a tangency point of $\mathcal{S}_{\sigma_0}$ with $\mathcal{C}_{\sigma_0,\, \Escr}\,$ we necessarily have that $\sigma_2=0$, as well as the tangency condition

\mathrm{rank}\begin{pmatrix}
\displaystyle{2\sigma_1} & \displaystyle{2\sigma_2} & \displaystyle{2\sigma_3} \\
\displaystyle{\frac{\partial \Hscr_{int}}{\partial \sigma_1}} & \displaystyle{0} & \displaystyle{\frac{\partial \Hscr_{int}}{\partial \sigma_3}}
\end{pmatrix}=1~~.

The latter condition implies that

\dot{\sigma}_2=\sigma_3\frac{\partial \Hscr_{int}}{\partial \sigma_1}-\sigma_1\frac{\partial \Hscr_{int}}{\partial \sigma_3}=0

that is, the point of tangency is a fixed point of the flow. Up to terms of second order in the variables $\sigma_i$ (i.e. of fourth order in the eccentricities), we find

\Hscr_{int}=A\sigma_1^2+B\sigma_3^2+C\sigma_1\sigma_3+D(\sigma_0)\sigma_1+E(\sigma_0)\sigma_3+F(\sigma_0)+...

where $A,B,C$ are constants, while the functions $D(\sigma_0)$ and $E(\sigma_0)$ are linear in $\sigma_0$ and $F(\sigma_0)$ contains terms linear and quadratic in $\sigma_0$. The quadratic form $A\sigma_1^2+B\sigma_3^2+C\sigma_1\sigma_3$ yields hyperbolas, being $A,B,C$ such that $C^2>AB\,$. Thus, for any permissible value $\sigma_0$ (or, equivalently, of the integral $J\,$), the surface $\mathcal{C}_{\sigma_0,\, \Escr}$ intersects the plane $\sigma_2=0$ along hyperbola-like curves (Fig.<ref>). The two points of tangency occur at the values $\sigma_0^{(A)}=2J^{(A)}$ and $\sigma_0^{(B)}=2J^{(B)}$. We find that $\sigma_0^{(A)}<\sigma_0^{(B)}$, while, checking the sign of $\cos(w_3-w_2)$ for the corresponding fixed points, we identify the left tangency (see Fig.<ref>) as the B-mode and the right tangency as the A-mode at the given level of energy.

With the help of the bottom-right panel of Fig.<ref> it is possible, now, to interpret the form of the phase portraits as in the left column of the same figure. To this end, we specify the correspondence between the various curves of the phase flow on the sphere, obtained by the intersections between the surfaces $\mathcal{S}_{\sigma_0}$ and $\mathcal{C}_{\sigma_0,\Escr}$ as $\sigma_0$ is altered in the interval $\sigma_0^{(A)}\leq\sigma_0\leq\sigma_0^{(B)}$, and the mapping of these curves to the surface of section $(\e_2\cos\omega_2,\e_2\sin\omega_2)$. Since $\Hscr_{int}$ does not depend on $\sigma_2$, all the curves produced by intersections of the surfaces $\mathcal{S}_{\sigma_0}$ and $\Cscr_{\sigma_0,\,\Escr}$ contain points which lie in the meridian circle produced by the intersection of $\mathcal{S}_{\sigma_0}$ with the plane $\sigma_2=0$. In particular, the points of tangency $A$ and $B$ belong to this meridian. Besides these points, there are two critical curves which separate domains of libration of the angle $\omega_2-\omega_3$ around the value 0 (mode A), or $\pi$ (mode B), from domains where the angle $\omega_2-\omega_3$ circulates.

By varying, now, the value of $\sigma_0$ in the interval $\sigma_0^{(A)}\leq\sigma_0\leq\sigma_0^{(B)}$ we progressively obtain curves on the sphere which pass from a librating domain around the fixed point $A$ to a circulating domain, and then to a librating domain around the fixed point $B$. The first such transition occurs at a value $\sigma_0^{(S)}$ where the curve corresponding to the intersection between $\mathcal{S}_{\sigma_0^{(S)}}$ and $\mathcal{C}_{\sigma_0^{(S)},\,\Escr}$ passes from the south pole $S$ of the sphere $\mathcal{S}_{\sigma_0^{(S)}}$. The coordinates of the south pole are $\sigma_1^{(S)}=\sigma_2^{(S)}=0$, $\sigma_3^{(S)}=-\sigma_0^{(S)}$, implying $J=-\Gamma$ or $W_2+W_3=-(W_2-W_3)$, hence $W_2=0$. This means a curve in the Poincaré section (thick purple) which crosses the origin $\e_2=0$. As clear from Fig.<ref>, we stress the well known fact that this curve means no real separatrix in the surface of section, generated by any kind of unstable periodic orbit, but it merely reflects the singularity induced by projecting an (all continuous) transition taking place in the phase space of the integrable Hamiltonian $\Hscr_{int}$, which is the 3-sphere, to the usual Poincaré section applicable to the full problem, i.e., the plane $(\e_2\cos\omega_2,\e_2\sin\omega_2)$.

Passing, now, the value $\sigma_0=\sigma_0^{(S)}$, we have curves of the sphere which are projected to invariant curves still surrounding the fixed point $A$, but for which the argument $\omega_2-\omega_3$ (or $\omega_2-\pi$, in the surface of section) circulates. A second limit of the circulation domain occurs at a value $\sigma_0^{(N)}$ where the curve corresponding to the intersection between $\mathcal{S}_{\sigma_0^{(N)}}$ and $\mathcal{C}_{\sigma_0{(N)}\, , \Escr}$ passes from the north pole $N$ of the sphere $\mathcal{S}_{\sigma_0^{(N)}}$. We readily find that the whole curve in the sphere $\mathcal{S}_{\sigma_0^{(N)}}$, except for the pole $N$ itself, yields the open black curve in the surface of section corresponding to the positive semi-plane $\e_2\cos(\omega_2)\geq 0$ (or $-\pi/2\leq \omega_2\leq \pi/2$), while the north pole itself inflates to the dotted semicircle obtained in the negative semi-plane $\e_2\cos(\omega_2)< 0$ (or $\pi/2< \omega_2<3\pi/2$). Here again the singularity is not real but only due to the choice of the variables representing the surface of section (see e.g. [Pauwels, 1983], [Henrard and Libert, 2004], [Michtchenko et al., 2006]). In fact, $N$ has coordinates $\left(\sigma_0^{(N)},\, \sigma_1^{(N)}=0,\,\sigma_2^{(N)}=0,\,\sigma_3^{(N)}=\sigma_0^{(N)}\right)\,$, implying $W_3=0$, i.e. $\e_3=0$. However, the equality $\sigma_3^{(N)}=\sigma_0^{(N)}$ implies also $X_2^2+Y_2^2=2\sigma_0$, i.e. a circle on the section $Y_3=0$. Together with the condition $\dot{Y}_3\geq 0$ of the Poincaré section, this implies the semi-circle $X_2\geq 0$, i.e. the dotted part of the black curve in the top-left panel of Fig.<ref>. It can be shown that the two parts or the curve join each other smoothly at two limiting values $X_2=0,Y_2=\pm Y_{2,max}=\pm \sqrt{2\sigma_0^{(N)}}$. In fact, the semi-circle $Y_2=\pm \sqrt{2\sigma_0^{(N)}-X_2^2}$, $X_2>0\,$, corresponding to the inflation of the north pole, yields $\lim_{X_2\rightarrow 0^+}(d Y_2/ d X_2)=0$, while the open curve $Y_2=Y_2(X_2)$, corresponding to all other points of the intersection of  $\mathcal{S}_{\sigma_0^{(N)}}$ and $\mathcal{C}_{\sigma_0{(N)}, \, \Escr}$ except for the pole, yields the limit

\begin{equation}
\label{derivata.curva.critica}
\lim_{X_2\rightarrow 0^-}\frac{d Y_2}{d X_2}=
\lim_{X_2\rightarrow 0^-}
\left(-\frac{\frac{\partial \Hscr_{int}}{\partial X_2}}
{\frac{\partial \Hscr_{int}}{\partial Y_2}}\right)~~.
\end{equation}

From the form of $H_{int}\,$, recalling Eq. (<ref>), we readily find $\partial\Hscr_{int}/\partial X_2=0$, $\partial\Hscr_{int}/\partial Y_2\neq 0$
at the north pole limit $X_2=X_3=Y_3=0$, $Y_2=\pm Y_{2,max}$.

Finally, for $\sigma_0$ in the interval $\sigma_0^{(N)}<\sigma_0\leq\sigma_0^{(B)}$ we find invariant curves in the sphere $\mathcal{S}_{\sigma_0}$ mapped to closed invariant curves around the tangency corresponding to the B-mode fixed point, which yield also closed curves in the Poincaré section for which the argument $\omega_2-\omega_3$ (or $\omega_2-\pi$ in the section) librates around the value $\pi$, i.e., $\omega_2$ librates around the value $\omega_2=0$ of the B-mode.

§.§.§ Semi-analytical (normal form) determination of the periodic orbits A and B

In the previous subsection we have seen how the existence of the A and B modes, which generalize the apsidal corotation periodic orbits of the model $\Hscr_{planar}$ to the spatial case, can be established within the framework of the integrable model $\Hscr_{int}=\Hscr_{planar}+\Hscr_{0,space}$. In the present subsection we discuss how to recover semi-analytically the periodic orbits A and B under the full Hamiltonian $\Hscr_{sec}=\Hscr_{int}+\Hscr_{1,space}$. Besides its relevance in establishing the existence of these orbits in the full model (up to an exponentially small remainder), a computation of the orbits A and B using normal forms allows to obtain a semi-analitycal representation of the long term time series of the orbital elements for these planetary trajectories. Since the modes A and B are among the most probable expected endstates of the formation process for exoplanetary systems, such a representation can be of use also in the interpretation of the observational data regarding the planetary orbital configurations in such systems.

Let $(\psi_*,\Gamma_*, J_\ast)$ be a fixed point of the integrable Hamiltonian $\Hscr_{int}(\psi,\Gamma;J)$, corresponding to one of the modes A or B. Introduce the translation

\begin{equation}\label{traaps}
\psi=\psi_{\ast}+\delta\psi,~~\Gamma=\Gamma_{\ast}+\delta\Gamma,~~J=J_{\ast}+\delta J.
\end{equation}

The transformation $(\psi,\varphi,\Gamma,J)\rightarrow (\delta\psi,\varphi,\delta\Gamma,\delta J)$ is canonical. The Hamiltonian $\Hscr_{sec}$ in the new variables reads:

\begin{equation}\label{hamdelta}
\Hscr_{sec}(\delta\psi,\varphi,\delta\Gamma,\delta J)=\Hscr_{int}(\delta\psi,\delta\Gamma,\delta J)
+ \Hscr_{1,space}(\delta\psi,\varphi,\delta\Gamma,\delta J)\, .
\end{equation}

We then have the following

Proposition: There is a near-to-identity canonical transformation

\begin{equation}\label{phiper}
(\delta\psi,\varphi,\delta\Gamma,\delta J)=\Phi^{(r)}(\delta\t{\psi},\t{\varphi},\delta\t{\Gamma},\delta\t{J})~~,
\end{equation}

real-analytic in a suitably defined domain, and obtained by a composition of $r$ Lie canonical transformations, with $r$ a sufficiently high integer, such that the Hamiltonian in the variables $(\delta\t{\psi},\t{\varphi},\delta\t{\Gamma},\delta\t{J})$ takes the form:

\begin{equation}\label{hamnophi}
\begin{split}
\Hscr^{(r)}(\delta\t{\psi},\t{\varphi},\delta\t{\Gamma},\delta\t{J})=
\end{split}
\end{equation}

with $||R^{(r)}||<<||Z^{(r)}||$ for a suitably defined norm in the domain of the transformation. The superscript $(r)$ in the above expressions means the expression derived after substituting the original variables $(\delta\psi,\varphi,\delta\Gamma,\delta J)$ with the new canonical variables obtained by the transformation inverse to the $r-$step transformation (<ref>):

\begin{equation}\label{phiper}
\left(\Phi^{(r)}\right)^{-1}(\delta\psi,\varphi,\delta\Gamma,\delta J)~~.
\end{equation}

The functions $Z^{(r)}$ and $R^{(r)}$ are hereafter called the `normal form' and `remainder' respectively.

The procedure by which to obtain the transformation $\Phi^{(r)}$ as well as its inverse is summarized in Appendix B, while a detailed demonstration of the above proposition is given elsewhere (Mastroianni and Efthymiopoulos, in preparation). In summary, this is a classical Birkhoff-like normalization procedure by which the angle $\varphi$ is eliminated from the Hamiltonian via a sequence of Lie canonical transformations (see [Efthymiopoulos, 2012]).

The physical meaning of the above normalization procedure is the following: in the planar-like regime, we have $||\Hscr_{1,space}||<<\Hscr_{int}$. Then, through the normalization we find a new set of variables in which, except for a very small remainder,  the dynamics locally (around the equilibrium $(\psi_*,\Gamma_*,J_*)$) is given by the normal form $Z^{(r)}(\t{\psi},\t{\Gamma},\t{J})$, with

\begin{equation}\label{tildevar}
\t{\psi}=\psi_*+\delta\t{\psi},~~\t{\Gamma}=\Gamma_*+\delta\t{\Gamma},~~\t{J}=J_*+\delta\t{J}~~.
\end{equation}

The phase flow induced by the integrable Hamiltonian $Z^{(r)}$ has the same structure as the one of the Hamiltonian $\Hscr_{int}$ analyzed in the previous subsection, differing by it just in the fact that $Z^{(r)}$ contains terms arising from the normalization of $\Hscr_{1,space}$ up to a very small remainder. In particular, Hamilton's flow under the normal form  in the transformed variables admits a periodic orbit given by

\begin{eqnarray}\label{nfper}
\t{\psi}_*&=&\psi_*+\delta\t{\psi}_*~~, \nonumber\\
\t{\Gamma}_*&=&\Gamma_*+\delta\t{\Gamma}_*~~, \\
\t{J}_*&=&J_*+\delta\t{J}_*~~, \nonumber\\
\t{\varphi}(t)&=&\t{\varphi}(0)+\nu_*t~~,\nonumber
\end{eqnarray}

where $(\t{\psi}_*$, $\t{\Gamma}_*,\t{J}_*)$ are computed by the system of algebraic equations

\begin{equation}\label{compnfper}
\begin{cases}
\displaystyle{\dot{\t{\psi}}=\frac{\partial Z^{(r)}}{\partial\t{\Gamma}}=0}
\\[2ex]
\displaystyle{\dot{\t{\Gamma}}=-\frac{\partial Z^{(r)}}{\partial \t{\psi}}=0}\\[1.3ex]
\end{cases}
\end{equation}

and the frequency $\nu_*$ is given by

\begin{equation}\label{nustar}
\nu_*=\left({\partial Z^{(r)}
\over\partial\t{J}}\right)_{\t{\psi}=\t{\psi}_*,\t{\Gamma}=\t{\Gamma}_*,\t{J}=\t{J}_*}~~.
\end{equation}

First row: Poincaré surfaces of section for the values of the energy (from left to right), $\Escr=-5.72\cdot 10^{-5},\, -3.19\cdot 10^{-5},\, -2.53\cdot 10^{-5}\,$. Second and third rows: time series of the evolution of the eccentricities for both planets along the mode $A$ (anti-apsidal, second row) and $B$ (apsidal, third row). The curves in red and in magenta show the time series for the eccentricities $\e_2$ and $\e_3$ respectively as computed semi-analytically. The black and blue curves, instead, show the time series $\e_2$ and $\e_3$ as computed through a numerical evaluation of the periodic orbits.

Figure <ref> shows an example of the comparison between the semi-analytical computation of the periodic orbits A, B on the basis of the normal form flow of $Z^{(r)}$, and a full numerical computation of the same orbits. The semi-analytical computation of the periodic orbits proceeds by the following steps:

i) we use the tangency method of subsection <ref> to compute the fixed points of the integrable Hamiltonian $\Hscr_{int}$, first in the Hopf variables and then in the variables $(\psi$, $\Gamma$, $J)\,$, obtaining $\psi_*$, $\Gamma_*$, $J_*$.

ii) Using an appropriate expansion of the Hamiltonian, as well as the method of composition of Lie series (see Appendix B), we then obtain the transformation $\Phi^{(r)}$ and its inverse $\left(\Phi^{(r)}\right)^{-1}$, as well as the normal form $Z^{(r)}$ representation of the full Hamiltonian.

iii) Implementing a Newton method, we then compute the root $(\t{\psi}_*$, $\t{\Gamma}_*,\t{J}_*)$ of the system of algebraic equations (<ref>), as well as the frequency (<ref>). This yields the time evolution of all
four quantities $(\delta\t{\psi}=\t{\psi}_*-\psi_*,\delta\t{\Gamma}=\t{\Gamma}_*-\Gamma_*,\delta\t{J}=\t{J}_*-J_*)$ (fixed) and $\t{\varphi}(t)$ as in Eq.(<ref>).

iv) We finally obtain the semi-analytical approximation to the time flow of all four variables $\psi(t)=\psi_*+\delta\psi(t)$, $\varphi(t)$, $\Gamma(t)=\Gamma_*+\delta\Gamma(t)$, $J(t)=J_*+\delta J(t)$, through the normalizing transformation

\begin{equation}\label{nfperflow}
(\delta\psi(t),\varphi(t),\delta\Gamma(t),\delta J(t))
=\left(\Phi^{(r)}(\delta\t{\psi},\t{\varphi},\delta\t{\Gamma},\delta\t{J})\right)_{ \delta\t{\psi}=\delta\t{\psi}_*,\t{\varphi}=\t{\varphi}(t), \delta\t{\Gamma}=\delta\t{\Gamma}_*,\delta\t{J}=\delta\t{J}_*}~~.
\end{equation}

This can be further transformed into a time series of the evolution of the orbital elements along the periodic orbit through the equations

\begin{align}\label{nfperele}
&\omega_2=\frac{-\psi-\varphi}{2}\, , & &W_2=J+\Gamma\, , \notag\\
& \omega_3=\frac{\psi-\varphi}{2}\, , & &W_3=J-\Gamma\,,
\end{align}

which allow to recover the evolution of the arguments of perihelia and eccentricities

\e_j=\sqrt{1-(1-W_j/\Lambda_j)^2},~~j=2,3

as well as (equivalently to Eq. (<ref>))

\begin{align}\label{nfpertheta}
\Theta_2&=\frac{(\AMD-W_2-W_3)(2 \Lambda_3+W_2-W_3-\AMD)}{2 (\Lambda_2+\Lambda_3-\AMD)}\\
\Theta_3&=\frac{(\AMD-W_2-W_3)(2 \Lambda_2-W_2+W_3-\AMD)}{2 (\Lambda_2+\Lambda_3-\AMD)}\, \notag
\end{align}

which allow to compute the time series for the inclinations

\i_j=\cos^{-1}\left(1-\frac{\Theta_j}{\Lambda_j-W_j}\right)

of the planetary periodic orbits.

As shown in Fig.<ref>, the semi-analytical representation of the periodic orbits fits with precision the numerical ones for both the modes A and B. As regards the physical interpretation, we note that the main effect of the perturbation $\Hscr_{1,space}$ is to induce a periodic oscillation in the eccentricities (and, hence, also the inclinations) of both planets, which are no longer constant, contrary to what holds for the classical apsidal corotation orbits in the planar case, or for the orbits of the integrable approximation $\Hscr_{int}$. From Fig.<ref> it is evident that the amplitude of the oscillation of the eccentricities of both planets increases with the energy $\Escr$, and, hence, with the mutual inclination. Most, notably, however, we point out the ability of the semi-analytical theory to well represent the orbits of both modes A and B in the regime after the onset of the saddle-node bifurcation giving rise to the periodic orbits $C_1$ and $C_2$, as in the third panel of Fig.<ref> (see subsection <ref>). This implies that, while at the corresponding level of energies the structure of the phase portraits is altered considerably with respect to the one of $\Hscr_{int}$, the representation of the motions by a local normal form computed in the neighborhood of each of the modes A or B yields the correct picture of the dynamics, both qualitatively and quantitatively.

§.§ Transition regime between planar-like and Lidov-Kozai:
sequences of bifurcations

The planar-like regime discussed in the previous subsection characterizes the structure of the phase portraits of the full system up to the energy $\Escr_C$, where the orbits $C_1,C_2$ are generated by a saddle node bifurcation in the taking place in the neighborhood of the mode B, as in the transition seen in the third row of Fig.<ref>. We call transition regime the one holding at energies in the interval $\Escr_{C}\leq\Escr_{C,2}$.

Poincaré surfaces of section for some values of the energy illustrating the sequence of bifurcations taking place in the `transition regime' (see text): $\Escr=-2.53\cdot 10^{-5}$  (top left), $\Escr=-2.08 \cdot 10^{-5}$ (top right), $\Escr=-1.9\cdot 10^{-5}$ (bottom left), $\Escr=-1.58\cdot 10^{-5}$ (bottom right), corresponding, respectively, to panels $11\,$, $13\,$, $14$ and $15$ of Figure <ref>.

Top row: the evolution of the orbital eccentricities $\e_2,\,\e_3$ for the periodic orbits corresponding to the modes A (top left panel: top curve $\e_2$, bottom curve $\e_3$), and B (top middle panel: top curve $\e_2$, bottom curve $\e_3$). The top right panel shows the librations of the argument $\omega_2-\omega_3$ for the B-mode (top curve) or the A-mode (bottom curve). Bottom row: the evolution of the orbital eccentricities $\e_2,\,\e_3$ for the periodic orbits corresponding to the fixed points $C_1$ (bottom left panel: top curve $\e_3$, bottom curve $\e_2$), and $C_2$ (bottom middle panel: top curve $\e_2$, bottom curve $\e_3$).The bottom right panel shows the librations of the argument $\omega_2-\omega_3$ for the orbit $C_2$ (top curve) or $C_1$ (bottom curve). The above time series are computed at the value of the energy $\Escr=-2.08\cdot 10^{-5}\,$, corresponding to the $13$-th panel of Figure <ref>.

Figure <ref> presents in detail the sequence of bifurcations taking place across the transition regime, which eventually lead to turning unstable the periodic orbit of mode B. At the beginning of the transition, for energies slightly larger than $\Escr_C$, the orbit $C_1$ turns from unstable to stable by a pitchfork bifurcation, which gives rise to an unstable periodic orbit yielding two fixed points, $D_1,D_2$, in the Poincaré surface of section. As the energy increases, these fixed points move initially away from the fixed point of the orbit $C_1$, while later (for still larger energy) they approach the fixed point of the mode B. Finally, at a second critical value of the energy the fixed points $D_1,D_2$ collide with the B-mode. This terminates the D-family of periodic orbits, by an inverse pitchfork bifurcation which renders the B-mode unstable.

Figure <ref> illustrates the evolution of the eccentricity vectors for all four periodic orbits $A,B,C_1,C_2$.
In all four cases the eccentricities of both planets oscillate periodically, while the argument $\omega_2-\omega_3$ undergoes small librations around one of the values $0$ or $\pi$.

§.§ Lidov-Kozai regime

We finally discuss the transition seen in the last row of Figure <ref>, in which the periodic orbit $C_2$ turns from stable to unstable via the Lidov-Kozai mechanism. This is accompanied by a large volume of trajectories around $C_2$ becoming chaotic. In the case of a test inner particle ($m_2=0$) at circular orbit $\e_2=0$, the mutual inclination is a preserved quantity, equal to the inclination of the test particle $\i_{mut}=\i_2$ (since $m_2=0$ the Laplace plane coincides with the constant orbital plane of the outer particle). Furthermore, the stability character of the so-called Lidov-Kozai state ( [Lidov, 1962], [Kozai, 1962], [Lidov and Ziglin, 1976],
[Ito and Ohtsuka, 2019]) depends only on the value of the inclination. In particular, in the quadrupolar approximation, the transition from stability to instability occurs at a critical inclination equal to $\cos^{-1}\sqrt{\frac{3}{5}}\sim 39^{\circ}.2\,$ (see [Naoz, 2016] for a review).

Here, instead, we will use a criterion analogous to the one of the classical Lidov-Kozai mechanism in order to obtain an estimate of the critical energy $\Escr_{C,2}$ at which the orbit $C_2$ turns from stable to unstable in the framework of the quadrupolar approximation to the secular Hamiltonian $\Hscr_{sec}$ in the full three body problem, i.e., for $m_2\neq 0$. The quadrupolar Hamiltonian, apart from constants, reads

\begin{equation}
\label{Ham.quadrupolar.c}
\Hscr_{quad}=-\frac{3\,\cgrav\,m_2\,m_3}{8\,a_3}\left(\frac{a_2}{a_3}\right)^2\frac{1}{\left(1-\e_3^2\right)^{3/2}}\Fscr_{quad}\, ,
\end{equation}

\begin{equation}
\label{Ham.quadrupolar}
\Fscr_{quad}=-\frac{1}{3}-\frac{\e_2^2}{2}+\frac{3}{2}\,\e_2^2\,\theta^2+\theta^2+\frac{5}{2}\,\e_2^2\left(1-\theta^2\right)\cos(2\,\omega_2)\, ,
\end{equation}

with $\theta=\cos(\i_2+\i_3)\,$. Following the process of Jacobi reduction, without any book-keeping control, ammounts to expanding the cosine of the mutual inclinations as $\cos(\i_2+\i_3)=\cos(\i_2)\cos(\i_3)-\sin(\i_2)\sin(\i_3)\,$, we can use expression (<ref>), that automatically cancel the dependence of the Hamiltonian on the mutual inclination. The Hamiltonian takes now the form (apart from constants):

\begin{equation}\label{hamquad}
\begin{split}
\Hscr_{quad}=&-\frac{3 \,a_2 \mathcal{G} \left(3 \,\e_2^2+2\right) \,m_3^3}{64 \,a_3^2 \left(1-\,\e_2^2\right) \sqrt{1-\,\e_3^2} \,m_2}-\frac{\,a_2^2 \mathcal{G} \left(3 \,\e_2^2+2\right) \,m_2 \,m_3}{32 \,a_3^3 \left(1-\,\e_3^2\right)^{3/2}}+\frac{3 \,a_2^3 \mathcal{G} \left(3 \,\e_2^4-\,\e_2^2-2\right) \,m_2^3}{64 \,a_3^4 \left(1-\,\e_3^2\right)^{5/2} \,m_3}\\
&+L_z^2\left(\frac{3 \,a_2 \left(3 \,\e_2^2+2\right) \,m_3}{32 \,a_3^3 \left(1-\,\e_2^2\right) \left(1-\,\e_3^2\right)^{3/2} \,m_0 \,m_2}+\frac{3 \,a_2^2 \left(3 \,\e_2^2+2\right) \,m_2}{32 \,a_3^4 \left(1-\,\e_3^2\right)^{5/2} \,m_0 \,m_3}\right)\\
&-\frac{3 \,a_2 \left(3 \,\e_2^2+2\right)\,L_z^4}{64 \,a_3^4 \mathcal{G} \left(1-\,\e_2^2\right) \left(1-\,\e_3^2\right)^{5/2} \,m_0^2 \,m_2 \,m_3}\\
&+\cos(2\,\omega_2)\Bigg[\frac{15 \,a_2 \mathcal{G} \,\e_2^2 \,m_3^3}{64 \,a_3^2 \left(1-\,\e_2^2\right) \sqrt{1-\,\e_3^2} \,m_2}-\frac{15 \,a_2^2 \mathcal{G} \,\e_2^2 \,m_2 \,m_3}{32 \,a_3^3 \left(1-\,\e_3^2\right)^{3/2}}+\frac{15 \,a_2^3 \mathcal{G} \,\e_2^2 \left(1-\e_2^2\right) \,m_2^3}{64 \,a_3^4 \left(1-\,\e_3^2\right)^{5/2} \,m_3}\\
&\phantom{+\cos(2\,\omega_2)\Bigg[}+L_z^2\left(-\frac{15 \,a_2 \,\e_2^2 \,m_3}{32 \,a_3^3 \left(1-\,\e_2^2\right) \left(1-\,\e_3^2\right)^{3/2} \,m_0 \,m_2}-\frac{15 \,a_2^2 \,\e_2^2 \,m_2}{32 \,a_3^4 \left(1-\,\e_3^2\right)^{5/2} \,m_0 \,m_3}\right)\\
&\phantom{+\cos(2\,\omega_2)\Bigg[}+\frac{15 \,a_2 \,\e_2^2\,L_z^4}{64 \,a_3^4 \mathcal{G} \left(1-\,\e_2^2\right) \left(1-\,\e_3^2\right)^{5/2} \,m_0^2 \,m_2 \,m_3}\Bigg]\, .
\end{split}
\end{equation}

It can be observed that the secular Hamiltonian at the quadrupolar level does not depend on the argument of the pericenter of the outer planet $\omega_3\,$, therefore the system is integrable (a fact known as the “happy coincidence”, see [Lidov and Ziglin, 1976] or [Naoz, 2016] for a review). In particular, the non-dependence of the Hamiltonian on $\omega_3$ implies that the eccentricity of the outer planet $\e_3$ is a conserved quantity.

Using $\e_3$ as a parameter, the Hamiltonian $\Hscr_{quad}$ can now be regarded as a one degree of freedom system. This can obtain a polynomial form by passing to the Poincaré variables $(X_2, Y_2)$ described by (<ref>). After such a substitution, the quadratic part of the above Hamiltonian is given by

\begin{equation}
\label{Ham.quad.quadrupolar}
\Hscr_{2,quad}(X_2, Y_2; G_3, L_z)=\frac{\mathfrak{a}(G_3, L_z)}{2}\, Y_2^2-\frac{\mathfrak{b}(G_3, L_z)}{2}\, X_2^2 \, ,
\end{equation}

where the coefficients $\mathfrak{b}$ and $\mathfrak{a}$ are functions of $G_3=L_3\sqrt{1-\e_3^2}=$const and of $L_z\,$; in particular

\begin{equation}
\label{ab}
\begin{split}
&\mathfrak{b}= \frac{3\, \mathcal{G}^2\, L_2^3 \,\left(3 \,G_3^2 - L_2^2 + L_z^2\right) \,m_0 \,m_3^7}{8 \,G_3^5\, L_3^3 \,m_2^3}\\
&\mathfrak{a}=-\frac{3\, \mathcal{G}^2\, L_2 \left(5\, G_3^4 - 4\, G_3^2\, L_2^2 + 3\, L_2^4 - 10 \,G_3^2\, L_z^2 -
8\, L_2^2\, L_z^2 + 5\, L_z^4\right) m_0\, m_3^7}{16\, G_3^5\, L_3^3\, m_2^3}\, .
\end{split}
\end{equation}

Now, in order to find the critical value of the energy $\Escr_{C,2}$, for fixed $L_z$, at which the periodic orbit $\e_3=const$, $\e_2=0$ becomes unstable,  it is sufficient to compute the eigenvalues of the Jacobian matrix of the Hamiltonian vector field:

\begin{equation}
\label{matrix.M}
\displaystyle{\frac{\partial \dot{X_2}}{\partial X_2}} & \displaystyle{\frac{\partial \dot{X_2}}{\partial Y_2}}\\[2.5ex]
\displaystyle{\frac{\partial \dot{Y_2}}{\partial X_2}}& \displaystyle{\frac{\partial \dot{Y_2}}{\partial Y_2}}
\end{pmatrix}=
\begin{pmatrix}
\displaystyle{\frac{ \partial^2 \Hscr_{quad}}{\partial X_2\,\partial Y_2}} & \displaystyle{\frac{\partial^2 \Hscr_{quad}}{\partial^2 Y_2}}\\[2.5ex]