\displaystyle{-\frac{\partial^2 \Hscr_{quad}}{\partial^2 X_2}}& \displaystyle{-\frac{\partial^2 \Hscr_{quad}}{\partial Y_2\,\partial X_2}}
\end{pmatrix}=
\begin{pmatrix}
\displaystyle{0} & \displaystyle{\mathfrak{a}}\\
\displaystyle{\mathfrak{b}}& \displaystyle{0}
\end{pmatrix}\, .
\end{equation}

The transition occurs at a critical value of $\e_3=\e_{3,C_2}$ at which the eigenvalues of $M$ pass from imaginary to real. The critical energy is then given by $\Escr_{C,2}=\Hscr_{quad}(\e_2=0,\e_{3,C_2};L_z)$. The critical value $\e_{3,C_2}$ can be computed by the following

Consider the secular Hamiltonian developed up to a quadrupolar expansion. Define the quantities

\begin{equation}
\begin{split}
\label{critical.points}
&A=\frac{1}{5}\left(4L_2^2+5 L_3^2(1-\e_3^2)-L_2\sqrt{L_2^2+60 L_3^2(1-\e_3^2)}\right)\, ,\\
&B=\frac{1}{5}\left(4L_2^2+5 L_3^2(1-\e_3^2)+L_2\sqrt{L_2^2+60 L_3^2(1-\e_3^2)}\right)\, ,\\
&C=L_2^2-3 L_3^2 (1-\e_3^2)\, ,
\end{split}
\end{equation}

with $L_2=m_2\sqrt{\mathcal{G} m_0 a_2}>0\,$, $L_3=m_3\sqrt{\mathcal{G} m_0 a_3}>0\,$.
Then, the following cases hold:

Case i): $0< L_2<\sqrt{3}L_3$ $\wedge$ $0\leq \e_3\leq \sqrt{1-\frac{L_2^2}{3 L_3^2}}\,$, the periodic orbit $C_2$ is Floquet-stable if $0<L_z^2<A\,$ or $L_z^2 > B\,$ and Floquet-unstable if $A<L_z^2<B\,$.

Case ii):
$0< L_2<\sqrt{3}L_3$ $\wedge$ $ \sqrt{1-\frac{L_2^2}{3 L_3^2}}< \e_3 < \sqrt{1-\frac{L_2^2}{4 L_3^2}}\,$ or $\sqrt{3} L_3 \leq L_2< 2 L_3$ $\wedge$ $0\leq \e_3 < \sqrt{1-\frac{L_2^2}{4 L_3^2}}\,$, the periodic orbit $C_2$ is Floquet-stable if $C<L_z^2<A\,$ or $L_z^2>B\,$  and Floquet-unstable if $0<L_z^2<C\,$ or $A<L_z^2<B\,$ .

Case iii):
$0< L_2\leq 2 L_3$ $\wedge$ $  \e_3 = \sqrt{1-\frac{L_2^2}{4 L_3^2}}\,$ , the periodic orbit $C_2$ is Floquet-stable if $L_z^2>B\,$ and Floquet-unstable if $0<L_z^2<B\,$ .

Case iv):
$0< L_2\leq 2 L_3$ $\wedge$ $ \e_3 > \sqrt{1-\frac{L_2^2}{4 L_3^2}}\,$ or if $ L_2> 2 L_3$ , the periodic orbit $C_2$ is Floquet-stable if $A<L_z^2<C\,$ or $L_z^2>B\,$ and Floquet-unstable if $0<L_z^2<A\,$ or $C<L_z^2<B\,$.

Moreover, having the critical points of $L_z^2\,$, it is easy to compute the critical values for the mutual inclination, being  (by Eq. (<ref>))

\begin{align}
\label{i.max}
&\cos(\i_2+\i_3)=\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}}\, , &
&\max \i_{mut}=\arccos\left(\frac{L_z^2-L_2^2-L_3^2}{2 \,L_2\, L_3}\right)\, .
\end{align}

The proof of the previous Proposition is reported in the Appendix <ref>. We readily find that in the limit $m_2\rightarrow 0$, $e_3\rightarrow 0$, the Kozai angles $\i=39.2^\circ $ and $\i=140.77^\circ$ are recovered (see [Naoz, 2016]). In fact, for $\e_3=0\,$, the critical values $A\,$, $B$ and $C$ become

\begin{align}
\label{val.crit.circolare}
&A=\frac{1}{5}\left(4L_2^2+5 L_3^2-L_2\sqrt{L_2^2+60 L_3^2}\right)\, , &
&B=\frac{1}{5}\left(4L_2^2+5 L_3^2+L_2\sqrt{L_2^2+60 L_3^2}\right)\, , & &C=L_2^2-3\,L_3^2\, .
\end{align}

Then, at the limit $m_2\rightarrow 0\,$ we readily obtain that

\begin{align*}
&\left(\cos i_{mut}\right)_{Lz^2=A}=\frac{A-L_2^2-L_3^2}{2 L_2 L_3}=-\frac{\sqrt{L_2^2+60\,L_3^2}+L_2}{10\,L_3}\;\; \xrightarrow{L_2\rightarrow 0}\;\; -\sqrt{\frac{3}{5}}\\
&\left(\cos i_{mut}\right)_{Lz^2=B}=\frac{B-L_2^2-L_3^2}{2 L_2 L_3}=\phantom{-}\frac{\sqrt{L_2^2+60\,L_3^2}-L_2}{10\,L_3} \;\; \xrightarrow{L_2\rightarrow 0}\;\; \phantom{-}\sqrt{\frac{3}{5}}\, .
\end{align*}

Comparison between the phase portraits (Poincaré surfaces of section) obtained with the quadrupolar approximation $\Hscr_{quad}$ (left column) and the octupolar approximation $\Hscr_{oct}$ (right column), at the energy levels $\Escr=-4.9\cdot 10^{-6}$ (top), and $\Escr=-6.9\cdot 10^{-7}$ (bottom). The quadrupolar approximation yields the phase portrait of an integrable system, which contains a figure-8 separatrix for energies beyond $\Escr_{C,2}$. The octupolar approximation is necessary to obtain both periodic orbits $C_1$ and $C_2$, as well as the chaotic zone around $C_2$ caused by the overlapping of resonances in the neighborhood of this orbit.

Note that, while the above proposition strictly establishes the limit of $\Escr_{C,2}$ only in the quadrupolar approximation, in practice we find that the estimate is quite precise when the full Hamiltonian with higher order multipoles is considered. For example, applying the above criterion in the numerical example of Fig.<ref> yields the estimate $\Escr_{C,2}\approx -2.2\cdot 10^{-6}$, which is in agreement with the value of the energy where the transition is observed in the numerical phase portraits. On the other hand, since the quadrupolar approximation yields an integrable Hamiltonian, for any value of the energy $\Escr>\Escr_{C,2}$ $\Hscr_{quad}$ yields a figure-8 phase portraits with invariant curves surrounding the `frozen' (stable) orbits on both sides of the unstable orbit $\e_2=0$. However, this picture changes by adding just the octupolar terms to the model, given by

\begin{equation}
\Hscr_{oct}=\Hscr_{quad}+\t\Hscr_{oct}\, ,
\end{equation}

\begin{equation}
\label{Ham.octupolar.c}
\t\Hscr_{oct}=\frac{75\,\cgrav\,m_2\,m_3}{64\,a_3}\left(\frac{a_2}{a_3}\right)^3\frac{1}{\left(1-\e_3^2\right)^{5/2}}\,\e_2\,\e_3\,\Fscr_{oct}
\end{equation}

\begin{equation}
\begin{split}
\label{Ham.octupolar}
\Fscr_{oct}\,=&\,\frac{1}{40}\left(3\,\e_2^2+4\right)\left(1+11\theta-5\theta^2-15\,\theta^3\right)\cos(\omega_2-\omega_3)\\
+& \frac{7}{8}\,\e_2^2(\theta^2-1)(1+\theta)\cos(3\omega_2-\omega_3)+\frac{7}{8}\,\e_2^2(\theta^2-1)(1-\theta)\cos(3\omega_2+\omega_3)\, ,
\end{split}
\end{equation}

with $\theta=\cos(\i_2+\i_3)\,$. Note that the above formula is equivalent, after some algebraic operations, to the formulas given in [Migaszewski and Goździewski, 2011] and [Naoz, 2016], see also [Lithwick and Naoz, 2011], [Naoz et al., 2013], [Ford et al., 2000]. As depicted in Fig.<ref>, the adition of more harmonics besides $\cos(2\omega_2)$ via the octupole approximation implies the creation of a large domain of resonance overlap, mostly between the islands around the frozen orbits and the island of the $C_1$ orbit, which appears already in the octupole approximation. This implies, in turn, the disappearance of most quasi-periodic trajectories around the frozen orbits, and the appearance, instead, of a large domain of chaotic orbits, in accordance to what is observed in the last two panels of Fig.<ref> for the full model.

Finally, we note that the use of the above proposition for the determination of critical values for the transition to instability has to be done in conjunction with the test that $\cos(\max i_{mut})$ does not obtain unphysical values $|\cos(i_{mut})|>1$. Such an unphysical value can occur in case ii), for the critical value $L_z^2=C\,$ and in case iv), for the critical value $L_z^2=A\,$. For example, in the second case of ii), setting $\sqrt{3}L_3\leq L_2< 2\, L_3$ and $\e_3=0\,$, we find $C=L_2^2-3\,L_3^2$ which, together with $L_z^2=C\,$, leads to $\cos(\max i_{mut})=-2\,L_3/L_2<-1$ iff $L_2< 2\,L_3$. Similarly, in the second case of iv), where $L_2>2\,L_3$, setting $\e_3=0\,$, by the relation $L_z^2=A=\frac{1}{5}\left(4L_2^2+5 L_3^2-L_2\sqrt{L_2^2+60 L_3^2}\right)\,$, we obtain $\cos(\max i_{mut})=-(\sqrt{L_2+60\,L_3^2}+L_2)/(10\,L_3)$, which is physical only iff $L_2\leq 2\,L_3\,$.

§.§ Where does the $\upsilon-$Andromedae system lie in phase space?

Poincaré surface of section for the energy $\Escr=-2.081\cdot 10^{-5}$. The red curve (1) shows the orbit obtained by adopting our basic assumption as regards the initial conditions of the $\upsilon-$Andromedae system (see subsection <ref>), i.e. (after reduction to the Laplace plane) $\e_2=0.2445$, $\e_3=0.316$, $\omega_2=289.049^\circ$, $\omega_3=235.464^\circ$. The black curve (2) shows the trajectory obtained by changing by only about 10% the eccentricity $\e_2$ and argument of the periastron $\omega_2$, while maintaining the energy and argument of periastron of the outer planet invariant ($\e_2=0.269$, $\omega_2=296.1^\circ$, $\omega_3=235.464^\circ$, $\e_3=0.299$).

We finally comment on the implications of the above analysis of the phase space structure on the interpretation of the data given by astronomical observations. Figure <ref> shows the phase portrait (surface of section) for the value of the energy $\Escr=-2.081\cdot 10^{-5}$, corresponding to the one found after reducing to the Laplace plane the data for the observed $\upsilon$-Andromedae system (see subsection <ref>). The invariant curve (1) (marked in red) corresponds to the actual initial conditions of the system. According to the figure, this is a quasi-periodic orbit surrounding the family $C_2$. We note, however, that the orbit is very close to the separatrix between the $C_2$ librational domain and the domain of stability surrounding the fixed point of the mode B orbit. As a consequence, a change of the initial conditions $(\e_2,\omega_2)$ by less than 10% results in a trajectory ((2), black curve) which undergoes quasi-periodic oscillations around the (aligned) apsidal corotation orbit of the system. The proximity of the real trajectory to one of the separatrices generated by the families $C_1,\,C_2$ is noteworthy, as it is in contrast with the basic scenario of evolution of planetary orbits, according to which the dissipative phase of planetary migration should lead to an endstate close to one of the stable fixed points of the system (see, for example, [Beaugé et al., 2012]).

§ PARAMETRIC STUDY

Poincaré surfaces of section in the plane $(\e_2 \cos(\omega_2 ), \e_2 \sin(\omega_2 ))$ with $L_z$ fixed such that max$\,i_{mut}$ is $45^\circ\,$ and
different values of energy. We consider the models $a_2/a_3=1/3$ and (from top to bottom) $m_2/m_3=1/10,\,1/3,\,1,\,3,\,10\,$. The values of the energy (from top to bottom) are (from left to right), top $\Escr=-2.71\cdot 10^{-6}, -1.01\cdot 10^{-6}, -9.24\cdot 10^{-7},-5.73\cdot 10^{-7}, -8.23\cdot 10^{-8} $, $\Escr=-9.03\cdot 10^{-7},-3.37\cdot 10^{-7},-2.87\cdot 10^{-7}, -1.76\cdot 10^{-7},  -2.35\cdot 10^{-8} $, $\Escr=-2.68\cdot 10^{-7},-1.39\cdot 10^{-7},-1.14\cdot 10^{-7},-6.37\cdot 10^{-8}, -6.98\cdot 10^{-9}  $, $\Escr=-1.12\cdot 10^{-6},-9.26\cdot 10^{-7},-3.37\cdot 10^{-7}, -2.08\cdot 10^{-7},-1.67\cdot 10^{-8}  $, bottom $\Escr= -2.75\cdot 10^{-6}, -1.91\cdot 10^{-6},-4.72\cdot 10^{-7},-3.54\cdot 10^{-7}, -2.91\cdot 10^{-8}\,$. The surfaces of section have been computed by a numerical integration of trajectories in a Hamiltonian averaged in closed form with a multipolar expansion truncated at degree $6$  and expanded up to order $12$ in the eccentricities.

The purpose of the present section is to provide an overview of how the form of the phase portraits (Poincaré sections) is varied with the energy $\Escr$, and hence with an increase in the level of mutual inclination, in various types of planetary systems differing from our so far main example as regards the choice of mass ratio $m_2/m_3$ as well as of semi-major axis ratio $a_2/a_3$ between the two planets. A thorough parametric study of the latter question is beyond our present scope. Instead, here we focus on only one central aspect of this study, namely the question of how generic is the description, as in section <ref>, of the transition from the planar-like to the Lidov Kozai regime, when either the mass ratio $m_2/m_3$ or the semi-major axis ratio $a_2/a_3$ are altered, including the so-called hierarchical limits, which correspond to mass ratio limits $m_2/m_3\rightarrow 0$ (restricted three-body problem with the inner planet being a test particle), or $m_2/m_3\rightarrow\infty$ (the outer planet is a test particle), and semi-major axis ratio $a_2/a_3<<1$ (outer planet way further from the star than inner planet). To produce a suite of numerical experiments covering most cases of practical interest, we consider in the sequel five different values of the mass ratio $m_2/m_3$ representative of all possible mass hierarchies (or lack thereof) namely $m_2/m_3=1/10,\,1/3,\,1,\,3,\,10$, produced by the corresponding combination of the masses $m_0=1M_\odot\,$ and $1\,M_{J}\,$, $3\,M_{J}$ or $10\,M_{J}$ for $m_2$ and $m_3\,$, while we consider the semimajor axes ratio $a_2/a_3=1/7$ ($a_2=1\,AU$, $a_3=7\,AU$) as a case representing distance hierarchy, or $a_2/a_3=1/3$ ($a_2=1\,AU$, $a_3=3\,AU$) as a case of no distance hierarchy. In each one of the previous cases, a comparison with the analysis of section <ref> requires keeping fixed the value of the $\AMD$ (or $L_z$) while computing phase portraits with a varying value of the energy $\Escr$. To find a relevant value for $L_z$ in each experiment, we consider a uniform value of the maximum possible mutual inclination $\i_{max}$ (Eq.<ref>) in all the experiments, set as $\i_{mut}=45^\circ$. This typically proves to be slightly above the limit of the Lidov-Kozai instability (which is equal to $39.2^\circ\,$ in the quadrupolar approximation of the restricted three body problem with $m_2=0$ and turns to be about $\sim 42.5^\circ\,$ in our main numerical example of section <ref>). By Eq.(<ref>), this implies setting

with $\Lambda_2,\Lambda_3$ specified by the choice of masses and semi-major axes in each experiment. Finally, in order to get some safety from any accuracy issues, we make computations with a Hamiltonian expanded at slightly larger orders with respect to those of the main numerical example of the previous section, i.e. multipole order 6, and order 12 in the eccentricities. To arrive at the final secular model, all Hamiltonians are processed, and their phase portraits computed, in the same way as in section <ref>.

Figure  <ref> shows the phase portraits obtained for increasing energy $\Escr$ in the case $a_2/a_3=1/3$, and in various mass hierarchy scenarios, namely, for the values $m_2/m_3$ (from top to bottom row) $1/10$, $1/3$,
$1/1$, $3/1$, $10/1$.  The values of the energies whose corresponding phase portraits are shown are selected in each panel of Fig.<ref> so as to visualize the most important changes observed in the phase-space structure as the energy increases.

Our main numerical example (Fig.<ref>) exhibits a similar structural change in the phase portraits, with increasing energy, as in the fully non-hierarchical case $m_2/m_3=1$ illustrated in the third row of Fig.<ref>. However, in Fig.<ref> it is clear that the transition from the planar-like regime, where the modes A and B dominate the phase space, through the saddle-node bifurcation giving rise to the orbits $C_1$ and $C_2$ (transition regime), and, finally, the Lidov-Kozai regime (orbit $C_2$ becomes unstable and surrounded by a chaotic figure-8 separatrix-like layer) is generic in all cases $m_2/m_3\leq 1$ (top three rows of Fig.<ref>) and follows by the same sequence of bifurcations. As the ratio $m_2/m_3$ tends to small values, the main differences observed with respect to the non-hierarchical case are: i) in the planar like regime (i.e. top row, first panel), the domain of quasi-periodic orbits around the aligned apsidal corotation orbit (B-mode) occupies most of the available phase space, and ii) at energies beyond the one of the Lidov-Kozai instability, there is a significant domain of regular orbits obtained around two stable periodic orbits of non-zero eccentricity $\e_2$ with fixed points along the axes $\omega_2=\pm\pi/2$ in the surface of section (the so-called `frozen orbits'). Both properties (i) and (ii) can be interpreted by analogy with the restricted three-body problem $(m_2=0)$, in which, for low inclinations, the only fixed point of the quadratic (Laplace-Lagrange) secular Hamiltonian is associated with an eccentric orbit of eccentricity $\e_2$ equal to the forced value induced by the outer perturber and pericenter aligned with the one of the outer perturber (see [Ferraz-Mello, 2007]). On the other hand, in the cases $m_2/m_3>1$ we do not observe at all the transition to the Lidov-Kozai regime. Again, this can be interpreted in analogy with the RTBP, in which the outer Lidov-Kozai instability ($m_3\rightarrow 0$) occurs at an inclination ($\approx 63^\circ$; see [von Zeipel, 1910]) much higher than the maximum mutual inclination $\i_{max}=45^\circ$ considered in our examples. The case of a very high mutual inclination has been considered in the literature (see for example [Hansen and Naoz, 2020] and references therein), but it appears of rather theoretical interest compared to available observations on the orbital inclinations in exoplanetary systems.

Poincaré surfaces of section in the plane $(\e_2 \cos(\omega_2 ), \e_2 \sin(\omega_2 ))$ with $L_z$ fixed such that max$\,i_{mut}$ is $45^\circ\,$ and different values of energy. We consider the models $a_2/a_3=1/7$ and (from top to bottom) $m_2/m_3=1/10,\,1/3,\,1,\,3,\,10\,$. The values of the energy (from top to bottom) are (from left to right), top $\Escr=-1.54\cdot 10^{-7},-6.71\cdot 10^{-8},-5.78\cdot 10^{-8},-2.46\cdot 10^{-8}, -6.15\cdot 10^{-9} $, $\Escr=-4.5\cdot 10^{-8},-2.05\cdot 10^{-8}, -1.7\cdot 10^{-8},-7.16\cdot 10^{-9}, -1.79\cdot 10^{-9} $, $\Escr=-2.75\cdot 10^{-8},-8.25\cdot 10^{-9},-7.01\cdot 10^{-9},-4.97\cdot 10^{-9},  -5.51\cdot 10^{-10} $, $\Escr=-7.5\cdot 10^{-8},-5.94\cdot 10^{-8},-3.69\cdot 10^{-8}, -2.5\cdot 10^{-8},  -1.38\cdot 10^{-9} $, bottom $\Escr=-2.22\cdot 10^{-7},-1.21\cdot 10^{-7},-8.96\cdot 10^{-8},-4.91\cdot 10^{-8}, -2.42\cdot 10^{-9}\,$. The surfaces of section have been computed by a numerical integration of trajectories in a Hamiltonian averaged in closed form with a multipolar expansion truncated at degree $6$  and expanded up to order $12$ in the eccentricities.

Fig.<ref>, now, shows how the phase portraits evolve with the energy when, in addition to altering the mass ratio $m_2/m_3$, the limit of distance hierarchy is also approached ($a_2/a_3=1/7$). As an overall observation, we note again that the full sequence of bifurcations leading from the planar-like to the Lidov-Kozai regime is realized only in the cases $m_2/m_3\leq 1$. However, the main difference, with respect to the case $a_2/a_3=1/3$, is that in the hierarchical in distance case the phase portraits contain many more regular orbits, with a cosiderably large domain of stability around the frozen orbits surviving the system's perturbations even at the mass ratio $m_2/m_3=1$ (third row of Fig.<ref>. This can be interpreted by the fact that the distance hierarchy brings the system closer to the dynamics of the integrable Hamiltonian $\Hscr_{quad}$ (Eq.(<ref>)), reducing the relative importance of perturbations including and beyond the octupolar one.

§ CONCLUDING REMARKS

As a main concluding remark on what was exposed in sections <ref> to <ref> above, we wish to point out again the twofold purpose of the present paper. On one hand, we aim to establish a convenient formalism for the study of the 3D planetary three body problem, allowing for a straightforward classification, visualization, and semi-analytical study of the most important phenomena appearing in the phase space of secular motions as the mutual inclination between the planetary orbits increases from zero to a value ($\sim 45^\circ$) beyond the onset of the Lidov-Kozai instability. On the other, we aim to provide ourselves a study of some of these phenomena in some detail. In particular, we focus on the connection between a chain of major transitions observed in the phase space structure as the mutual inclination increases, and the associated bifurcations (and/or the change of stability character) of the most important periodic orbits of the problem. In physical terms, the hereby proposed formalism serves to analyze the gradual transition of the system between two well distinct regimes. One of them is a nearly-integrable regime, called above `planar-like', dominated by periodic orbits (A and B modes) analogous to the apsidal corotations of the planar case, as well as their surrounding quasi-periodic orbits. The second is the regime called above `Lidov-Kozai'. This is characterized by the highly inclined periodic orbits emerging in the octupolar approximation (orbits $C_1$ and $C_2$), one of which ($C_2$) undergoes the well known Lidov-Kozai transition from stable to unstable beyond a critical value of the mutual inclination.

As also emphasized in the introduction, our presently proposed formalism is strongly motivated by the need to provide a unifying conceptual framework in which to understand and put into context several results, distinct or partly overlapping with those presently reported, obtained in past literature by use of quite diverse methods and representations of the problem under exam. The following is a summary of our main conclusions, serving also as a navigation guide to the main methodological steps proposed in the preceeding sections of the paper.

i) Section <ref> exposed the methodology for obtaining a suitable secular Hamiltonian model, as well as our choice of variables and representation for phase portraits, aiming to present in a unified way both the planar-like regime, where apsidal corotation orbits dominate, and the Lidov-Kozai regime, where nearly circular highly inclined orbits dominate. The use of a particular book-keeping method was explained in subsection <ref>, leading, after Jacobi reduction, to a Hamiltonian decomposition of the form $\Hscr_{sec}=\Hscr_{planar}+\Hscr_{space}$, where the term $H_{planar}$ is integrable and the term $\Hscr_{space}$ depends on the system's $\AMD$ (or, equivalently, the value of the system's angular momentum $L_z$). We argue on the merits stemming from a choice of simple definition for the Poincaré surface section as in <ref>, allowing for a straightforward mapping of phase portraits to the evolution of the corresponding orbits in terms of Keplerian elements. This is accompanied by a detailed analysis of the permissible domains of initial conditions on the surface of section as the energy $\Escr$ is altered. It is demonstrated that these limits are defined by the intersections between any two of three important manifolds, namely the manifold of constant energy $\mathcal{M}(\Escr)$, the manifold of zero mutual inclination ${\cal{I}}_{0}$, and the condition of tangency of an orbit with the surface of section ($\dot{Y}_3=0$ in the Poincaré variables (<ref>)). In particular, we explain how the tangencies between $\mathcal{M}(\Escr)$ and ${\cal{I}}_0$ at two particular values of the energy $\Escr_A\,(=\Escr_{min})$, $\Escr_B\,(=\Escr_{2,3})$ mark the birth of the planar anti-aligned and aligned apsidal corotation orbits respectively. Finally, we test the precision of various models defined either by a closed-form scissors-averaged model obtained from a multipole expansion of the Hamiltonian truncated at degree $N_P$, or by the usual (Jacobi-reduced) Laplace-Lagrange expansion. In both the previous cases, we introduce an additional `book-keeping' order of truncation $N_{bk}$, with a book-keeping collecting at the same time powers of the eccentricities, of the mutual inclination and of the system's $\AMD$. As discussed in subsection <ref>, we find that, in general, a multipole expansion with $N_P=5$ or $6$ is sufficient to represent even non-hierarchical systems with semi-major axis ratio between the two planets $a_2/a_3=1/3$ or smaller. The multipole model has the salient features of being easy to scissors-average in closed form, while its lowest order truncations (quadrupolar and octupolar) yield formulas sufficiently small to offer analytical insight into the results. Comparison with the Laplace-Lagrange model (harder to obtain) shows the practical equivalence of the phase portraits at $N_P=5$ and with order of expansion as high as $N_{bk}=10$. On the other hand, obtaining the canonical transformation corresponding to the averaging, or going to second order in the masses, seems quite harder to achieve with the closed-form approach.

ii) Section <ref> gives a qualitative overview of the main mechanism behind the transition from the `planar-like' to the `Lidov-Kozai' regime, i.e. the birth of the families $C_1,C_2$ via a saddle-node bifurcation around the aligned AC (mode B), accompanied by the transition of the family $C_1$ from unstable to stable via a pitchfork bifurcation, and, finally, the de-stabilization of the B-mode via an inverse pitchfork transformation. These phenomena are described qualitatively via the corresponding phase portraits, i.e., surfaces of section obtained numerically in a basic example by increasing the value of the energy $\Escr$. We demonstrate how this latter increase is connected to the overall rise of the level of mutual inclination of the orbits (Fig.<ref>). Finally, we numerically observe the Lidov-Kozai instability in the same phase portraits which turns the family $C_2$ from stable to unstable at a critical energy $\Escr_{C,2}$.

iii) Based on the form of the phase portraits, we identify three basic regimes in which a two-planet system with inclined orbits can be found:

- the `planar-like' regime is characterized by a phase-portrait similar to the one of an integrable model $\Hscr_{int}$ subsection <ref>, which is the 3D analogue of the integrable model $\Hscr_{planar}$. As pointed out in past literature (see references in subsection <ref>) this is a model whose natural phase space is the two-sphere, represented by the Hopf variables of Eq.(<ref>). We first implement a tangency method ( [Cushman and Bates, 1997], [Palacián et al., 2006], [Marchesiello and Pucacco, 2016]) to locate the periodic orbits of the integrable secular model. Then, by a normal form approach, we demonstrate how the perturbations to the integrable model can be normalized and absorbed in the definition of the periodic orbits corresponding to the apsidal corotations in the non-planar case. The main conclusion is that both the A (anti-aligned) and B (aligned) modes of the integrable model turn into orbits whose eccentricities are not exactly constant, but undergo periodic oscillations with a secular period easy to compute semi-analytically. At the same time, the pericenters of the planets in either mode do not remain locked exactly to the alignment or anti-alignment condition. Instead, the difference $\omega_2-\omega_3$ also undergoes oscillations around the central value (0 for anti-alignment, $\pi$ for alignment). The amplitude of these oscillations grows with the energy, and hence the level of mutual inclination,  where each of these periodic orbits is computed.

- the `transition regime', starting with the birth of the families $C_1$, $C_2$ and ending with the termination of the  intermediate ($D_{1,2}$) family at a collision with the B-mode rendering the latter unstable

- the`Lidov-Kozai' regime, in which the structure of the phase space is similar to the one provided by the octupolar model. As well known in literature, the turning of the orbit $C_2$ from stable to unstable can be predicted analytically in the framework of the (integrable) quadrupolar approximation. We adapt this prediction to a proposition (subsection <ref>) allowing to compute approximately the value of the transition energy $\Escr_{C,2}$ and of the corresponding critical inclination for a system with fixed $\AMD$. Actually, a careful investigation of the same proposition shows that it can be applied also independently of the choice of parameter, i.e., using the outer planet's eccentricity (which is constant in the context of the quadrupolar model) as well as altering the value of the systems $\AMD$ (or angular momentum). At any rate, the result obtained by such an analysis is only approximate, at least for non-hierarchical systems, since the octupolar terms introduce a significant perturbation to the quadrupolar model. It is noteworthy that the effect of the octupolar terms can be locally absorbed by a perturbation theory eliminating all harmonics except for $\cos(2\omega_2)$ (attempts to provide such a theory in the framework of the restricted three body problem are nearly as old as the subject itself, see [Ito and Ohtsuka, 2019] for a review). We are aware of no work attempting to provide the thresholds for the onset of the Lidov-Kozai instability in the framwework of the general three-body problem at truncations $N_P\geq 3$ by analytical means (normal forms), thus, this is proposed as a subject for further study.

iv) We give an overview of how the picture of the phase-space structure and dynamics depicted as above applies to systems in which we have either a mass hierarchy, or a distance hierarchy. As regards models with mass hierarchy, we give numerical evidence that a necessary condition to obtain the full chain of transitions from the `planar-like' to the `Lidov-Kozai' regime is that $m_2$ be smaller than $m_3$. This is justified by the fact that the critical mutual inclination for obtaining the Lidov-Kozai instability is $\approx 39^\circ$ in the limit $m_2<<m_3$, but it raises, instead, to $\approx 63^\circ$ in the limit $m_3<<m_2$. Such a high value of the mutual inclination is rather improbable in exoplanetary systems. On the other hand, a quick comparison of Figs.<ref> and <ref> shows that the so-called `planar-like' regime corresponds to a dynamical characterization (or, simply, the similarity of phase portraits with the planar case) rather than a really small value of the mutual inclination. For example, in our basic model the planar-like regime persists up to the energy $\Escr_C$ of birth of the $C_1,C_2$ families, which corresponds to trajectories whose mutual inclination can be as high as $35^\circ$. On the other hand, as can be deduced from Fig.<ref>, in models with distance hierarchy we arrive rather quickly to forms of the phase portraits reminiscent of those obtained in the Lidov-Kozai regime, i.e., with the $C_1$ and $C_2$ orbits rather, than the mode A and B orbits, dominating the dynamics. Insight to this phenomenon is provided by the fact that the octupolar model gives a far better approximation to the full Hamiltonian in distance-hierarchical cases. In that model, however, the harmonics $\cos(\omega_2-\omega_3)$, which is the basis for the appearance of apsidal corotation orbits, has smaller amplitude than the harmonics $\cos(2\omega_2)$ (already present at quadrupolar level). This implies that the dominant periodic orbits are $C_1$ and $C_2$, as well as the frozen orbits arising after the Lidov-Kozai instability. Note that the orbits $C_1$ and $C_2$ have aligned and anti-aligned apsides, while, by the position of the corrsponding fixed points in the phase portraits, we readily deduce that they have smaller eccentricities than the one of the A mode (which remains stable at all values of the energy). Thus, once again, the difference between these periodic orbits is dynamical rather than morphological, i.e., they play distinct roles in shaping the phase space locally around them.

Acknowledgements: C.E. acknowledges the support of MIUR-PRIN 20178CJA2B ‘New frontiers of Celestial Mechanics: theory and applications’.

[Beaugé et al., 2003]
C. Beaugé, S. Ferraz-Mello, and T. Michtchenko.
Extrasolar planets in mean-motion resonance: apses alignment and
asymmetric stationary solutions.
The Astrophysical Journal, 5930 (2):0
1124–1133, 2003.

[Beaugé et al., 2012]
C. Beaugé, S. Ferraz-Mello, and T. A. Michtchenko.
Multi-planet extrasolar systems-detection and dynamics.
Research in Astronomy and Astrophysics, 120
(8):0 1044, 2012.

[Brouwer, 1959]
D. Brouwer.
Solution of the problem of artificial satellite theory without drag.
The Astronomical Journal, 64:0 378–397, 1959.
ISSN 0004-6256.
URL <https://doi.org/10.1086/107958>.

[Cavallari and Efthymiopoulos, 2022]
I. Cavallari and C. Efthymiopoulos.
Closed-form perturbation theory in the restricted three-body problem
without relegation.
Celestial Mechanics and Dynamical Astronomy, 134, 2022.

[Cushman and Bates, 1997]
R. H. Cushman and L. M. Bates.
Global aspects of classical integrable systems.
Birkhauser, 1997.
ISBN 0817654852.

[Deitrick et al., 2015]
R. Deitrick, R. Barnes, B. McArthur, T. R. Quinn, R. Luger, A. Antonsen, and
G. F. Benedict.
The three-dimensional architecture of the $\upsilon$ Andromedae
planetary system.
The Astrophysical Journal, 7980 (1):0 46,

[Efthymiopoulos, 2012]
C. Efthymiopoulos.
Canonical perturbation theory; stability and diffusion in
Hamiltonian systems: applications in dynamical astronomy.
In Workshop series of the asociacion argentina de astronomia,
volume 3, pages 3–146, 2012.

[Ferraz-Mello, 2007]
S. Ferraz-Mello.
Canonical perturbation theories: degenerate systems and
resonance, volume 345.
Springer Science & Business Media, 2007.

[Ford et al., 2000]
E. B. Ford, B. Kozinsky, and F. A. Rasio.
Secular evolution of hierarchical triple star systems.
The Astrophysical Journal, 5350 (1):0
385–401, 2000.

[Gronchi and Milani, 1998]
G. F. Gronchi and A. Milani.
Averaging on Earth-Crossing Orbits.
Celestial Mechanics and Dynamical Astronomy, 710
(2):0 109–136, 1998.

[Hadjidemetriou, 1975]
J. D. Hadjidemetriou.
The continuation of periodic orbits from the restricted to the
general three-body problem.
Celestial mechanics, 120 (2):0 155–174, 1975.

[Hadjidemetriou, 2006]
J. D. Hadjidemetriou.
Symmetric and asymmetric librations in extrasolar planetary systems:
a global view.
Celestial Mechanics and Dynamical Astronomy, 950
(1):0 225–244, 2006.

[Hansen and Naoz, 2020]
B. M. Hansen and S. Naoz.
The stationary points of the hierarchical three-body problem.
Monthly Notices of the Royal Astronomical Society,
4990 (2):0 1682–1700, 2020.

[Henrard and Libert, 2004]
J. Henrard and A.-S. Libert.
The secular planetary three body problem revisited.
Proceedings of the International Astronomical Union,
20040 (IAUC197):0 49–54, 2004.

[Ito and Ohtsuka, 2019]
T. Ito and K. Ohtsuka.
The Lidov-Kozai Oscillation and Hugo von Zeipel.
Monographs on Environment, Earth and Planets, 7:0
1–113, 11 2019.

[Jacobi, 1842]
C. G. J. Jacobi.
Sur l'élimination des noeuds dans le problème des trois
corps. Par M. Jacobi.
Astronomische Nachrichten, 20:0 81–102, 1842.

[Kozai, 1962]
Y. Kozai.
Secular perturbations of asteroids with high inclination and
The Astronomical Journal, 670 (9):0 591–598,

[Laskar and Robutel, 1995]
J. Laskar and P. Robutel.
Stability of the planetary three-body problem.
Celestial Mechanics and Dynamical Astronomy, 620
(3):0 193–217, 1995.

[Laughlin et al., 2002]
G. Laughlin, J. Chambers, and D. Fischer.
A dynamical analysis of the 47 Ursae Majoris planetary system.
The Astrophysical Journal, 5790 (1):0
455–467, 2002.

[Lee and Peale, 2003]
M. H. Lee and S. Peale.
Secular evolution of hierarchical planetary systems.
The Astrophysical Journal, 5920 (2):0
1201–1216, 2003.

[Libert and Delsate, 2012]
A.-S. Libert and N. Delsate.
Interesting dynamics at high mutual inclination in the framework of
the Kozai problem with an eccentric perturber.
Monthly Notices of the Royal Astronomical Society,
4220 (3):0 2725–2736, 2012.

[Libert and Henrard, 2007]
A.-S. Libert and J. Henrard.
Exoplanetary systems: The role of an equilibrium at high mutual
inclination in shaping the global behavior of the 3-D secular planetary
three-body problem.
Icarus, 1910 (2):0 469–485, 2007.
ISSN 0019-1035.

[Libert and Tsiganis, 2009]
A.-S. Libert and K. Tsiganis.
Kozai resonance in extrasolar systems.
Astronomy & Astrophysics, 4930 (2):0
677–686, 2009.

[Libert et al., 2011]
A.-S. Libert, C. Hubaux, and T. Carletti.
The global symplectic integrator: an efficient tool for stability
studies of dynamical systems. Application to the Kozai resonance in the
restricted three-body problem.
Monthly Notices of the Royal Astronomical Society,
4140 (1):0 659–667, 2011.

[Lidov, 1962]
M. Lidov.
The evolution of orbits of artificial satellites of planets under the
action of gravitational perturbations of external bodies.
Planetary and Space Science, 90 (10):0
719–759, 1962.

[Lidov and Ziglin, 1976]
M. L. Lidov and S. Ziglin.
Non-restricted double-averaged three body problem in Hill's case.
Celestial Mechanics, 130 (4):0 471–489, 1976.

[Lithwick and Naoz, 2011]
Y. Lithwick and S. Naoz.
The eccentric Kozai mechanism for a test particle.
The Astrophysical Journal, 7420 (2):0 94,

[Locatelli and Giorgilli, 2000]
U. Locatelli and A. Giorgilli.
Invariant tori in the secular motions of the three-body planetary
Celestial Mechanics and Dynamical Astronomy, 780
(1):0 47–74, 2000.

[Marchesiello and Pucacco, 2016]
A. Marchesiello and G. Pucacco.
Bifurcation sequences in the symmetric 1:1 Hamiltonian resonance.
International Journal of Bifurcation and Chaos, 260
(04):0 1630011, 2016.

[McArthur et al., 2010]
B. E. McArthur, G. F. Benedict, R. Barnes, E. Martioli, S. Korzennik, E. Nelan,
and R. P. Butler.
New observational constraints on the $\upsilon$ Andromedae system
with data from the Hubble Space Telescope and Hobby-Eberly
The Astrophysical Journal, 7150 (2):0
1203–1220, 2010.

[Michtchenko et al., 2006]
T. Michtchenko, S. Ferraz-Mello, and C. Beaugé.
Modeling the 3-D secular planetary three-body problem: Discussion
on the outer $\upsilon$ Andromedae planetary system.
Icarus, 1810 (2):0 555–571, 2006.

[Michtchenko and Malhotra, 2004]
T. A. Michtchenko and R. Malhotra.
Secular dynamics of the three-body problem: application to the
$\upsilon$ Andromedae planetary system.
Icarus, 1680 (2):0 237–248, 2004.

[Migaszewski and Goździewski, 2008]
C. Migaszewski and K. Goździewski.
A secular theory of coplanar, non-resonant planetary system.
Monthly Notices of the Royal Astronomical Society,
3880 (2):0 789–802, 2008.

[Migaszewski and Goździewski, 2011]
C. Migaszewski and K. Goździewski.
The non-resonant, relativistic dynamics of circumbinary planets.
Monthly Notices of the Royal Astronomical Society,
4110 (1):0 565–583, 2011.

[Murray and Dermott, 1999]
C. D. Murray and S. F. Dermott.
Solar system dynamics.
Cambridge university press, 1999.

[Naoz, 2016]
S. Naoz.
The eccentric Kozai-Lidov effect and its applications.
Annual Review of Astronomy and Astrophysics, 54:0
441–489, 2016.

[Naoz et al., 2013]
S. Naoz, W. M. Farr, Y. Lithwick, F. A. Rasio, and J. Teyssandier.
Secular dynamics in hierarchical three-body systems.
Monthly Notices of the Royal Astronomical Society,
4310 (3):0 2155–2171, 2013.

[Palacián et al., 2006]
J. F. Palacián, P. Yanguas, S. Fernández, and M. A. Nicotra.
Searching for periodic orbits of the spatial elliptic restricted
three-body problem by double averaging.
Physica D: Nonlinear Phenomena, 2130 (1):0
15–24, 2006.

[Pauwels, 1983]
T. Pauwels.
Secular orbit-orbit resonance between two satellites with non-zero
Celestial mechanics, 300 (3):0 229–247, 1983.

[Petit et al., 2018]
A. C. Petit, J. Laskar, and G. Boué.
Hill stability in the AMD framework.
Astronomy & Astrophysics, 617:0 A93, 2018.

[Szeto and Lambeck, 1982]
A. Szeto and K. Lambeck.
On eccentricity functions for eccentric orbits.
Celestial mechanics, 270 (4):0 325–337, 1982.

[Thomas and Morbidelli, 1996]
F. Thomas and A. Morbidelli.
The Kozai resonance in the outer solar system and the dynamics of
long-period comets.
Celestial Mechanics and Dynamical Astronomy, 640
(3):0 209–229, 1996.

[von Zeipel, 1910]
H. von Zeipel.
Sur l'application des séries de M. Lindstedt à
l'étude du mouvement des comètes périodiques.
Astronomische Nachrichten, 183:0 345, 1910.

§ LAPLACE COEFFICIENTS

In this section we show how to compute the Laplace coefficients (calculated in Subsection <ref> numerically through equations (<ref>)) via a multipolar expansion. We can state the following

The Laplace coefficients, allowing to compute the direct part of the Hamiltonian (<ref>), i.e. (see (<ref>))

$$ \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(\alpha\right)\cos(j(\lambda_2-\lambda_3))\, ,$$

can be computed as

\begin{equation}
\label{Lapl.coeff.analitic}
\begin{split}
b_{s+\frac{1}{2}}^{(0)}\left(\alpha\right)&= 1+\sum_{k=1}^{+\infty}(-1)^{k}\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\,\alpha^{2k}\\
&\phantom{=}+\sum_{h\in\Ascr_2}\sum_{k=h}^{\infty}(-1)^{k-h} \frac{(s+\frac{1}{2})(s+\frac{3}{2})\ldots(s-\frac{1}{2}+k)}{\left(\frac{h}{2}\right)!\left(\frac{h}{2}\right)!(k-h)!}\alpha^{2k-h}\, ,\\
b_{s+\frac{1}{2}}^{(j)}\left(\alpha\right)&=\sum_{h\in\Ascr_j}\sum_{k=h}^{\infty}2(-1)^{k-h} \frac{(s+\frac{1}{2})(s+\frac{3}{2})\ldots(s-\frac{1}{2}+k)}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!(k-h)!}\alpha^{2k-h}\, ,\quad j\geq 1\, ,
\end{split}
\end{equation}

where $\alpha=a_2/a_3\,$ and $\Ascr_j=\big\lbrace h\in\naturali\,:\,h=\begin{cases}
2i\phantom{+1\:\,}\quad\textit{if j is even}\\
2i+1\quad\textit{if j is odd}
\end{cases} ,\,i\in\naturali\,,\: i\geq \lfloor \frac{j}{2}\rfloor \bigg\rbrace\, $, $j\geq 1\,$.

We are interested in the expansion of the following quantity

\begin{equation}
\label{start}
\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^{-(2s+1)}\frac{1}{\left(1+\alpha^2-2\alpha \cos(\lambda_2 - \lambda_3)\right)^{\frac{2\,s+1}{2}}}\, ,
\end{equation}

with $\alpha=a_2/a_3\,$; for simplicity, we can define $l=2s+1$ and $\sigma=\lambda_2-\lambda_3\,$ and we can Taylor expand the previous quantity[It is sufficient to think $2\alpha\cos(\sigma)-\alpha^2:=\varepsilon$ and to Taylor expand, in $\varepsilon=0\,$, the function $(1-\varepsilon)^{-l/2}\,$.], having

\begin{equation}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}=1+\sum_{k=1}^{+\infty}\frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\,\left(2\alpha \cos(\sigma)-\alpha^2\right)^k\, .
\end{equation}

Moreover, we can expand the quantity $\left(2\alpha \cos(\sigma)-\alpha^2\right)^k\,$, obtaining:

\begin{equation}
\begin{split}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}=&1+\sum_{k=1}^{+\infty}\frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\,\left(2\alpha \cos(\sigma)-\alpha^2\right)^k\, \\
=& 1+\sum_{k=1}^{+\infty}\frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\,\left[\sum_{h=0}^{k}\frac{k!}{h!(k-h)!}(-1)^{k-h}\alpha^{2(k-h)}(2\alpha\cos(\sigma))^h\right]\,\\
=& 1+\sum_{k=1}^{+\infty}\sum_{h=0}^{k}(-1)^{k-h}\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{h!(k-h)!}\alpha^{2k-h}\,2^h\,\cos(\sigma)^h\,.
\end{split}
\end{equation}

After having explicitely written the sum corresponding to $h=0\,$, it is possible to reverse the order of the sums (over $k$ and $h$), arriving to

\begin{equation}
\label{intermedia}
\begin{aligned}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}&=1+\sum_{k=1}^{+\infty}(-1)^{k}\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\alpha^{2k} \\
&\phantom{=}+\sum_{h=1}^{+\infty}\sum_{k=h}^{+\infty}(-1)^{k-h}\,2^h\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{h!(k-h)!}\alpha^{2k-h}\,\cos(\sigma)^h \\
&:= \sum_{h=0}^{+\infty}D_{h}\cos(\sigma)^h\, ,
\end{aligned}
\end{equation}

where we can define $D_h\,$ (with $h\geq 0$) as:

\begin{equation}
\label{Dh} D_h=
\begin{cases}
\displaystyle{1+\sum_{k=1}^{+\infty}(-1)^{k}\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{k!}\,\alpha^{2k}}\quad\quad\:\:\:\textit{for}\quad h=0\\
\displaystyle{\sum_{k=h}^{+\infty}(-1)^{k-h}\,2^h\, \frac{\frac{l}{2}(\frac{l}{2}+1)\ldots(\frac{l}{2}+k-1)}{h!(k-h)!}\,\alpha^{2k-h}}\quad\textit{for}\quad h\geq 1
\end{cases}.
\end{equation}

Therefore, in order to compute the Laplace coefficients, we need to expand also the $\cos(\sigma)^h\,$ in (<ref>); indeed, from one hand (Eq. (<ref>)) we have

$$ \frac{1}{\left(a_2^2 + a_3^2 - 2 \,a_2\,a_3\, \cos(\sigma)\right)^\frac{l}{2}}=a_3^{-l}\sum _{h\geq 0} D_h\cos(\sigma)^h\, ,

and on the other one (Eq. (<ref>))

$$ \frac{1}{\left(a_2^2 + a_3^2 - 2 \,a_2\,a_3\, \cos(\sigma)\right)^\frac{l}{2}}=a_3^{-l}\sum _{j\geq 0} b^{(j)}_\frac{l}{2}(\alpha)\cos(j\sigma)\, .

Then, using again the binomial formula, we have

\begin{equation}
\label{esp.coseno}
\cos(\sigma)^h=\left(\frac{e^{i\sigma}+e^{-i\sigma}}{2}\right)^h=\frac{1}{2^h}\sum_{k=0}^{h} \binom{h}{k}e^{ik\sigma} e^{-i(h-k)\sigma}=\frac{1}{2^h}\sum_{k=0}^{h} \binom{h}{k}\cos((2k-h)\sigma),
\end{equation}

where we applied the parity of the cosine. Then, inserting the previous expression (<ref>) in Eq. (<ref>), we obtain

\begin{equation}
\label{final0}
\begin{aligned}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}=&D_0+\sum_{h=1}^{+\infty}\sum_{k=0}^{h} \frac{D_{h}}{2^h}\frac{h!}{k!(h-k)!}\cos((h-2k)\sigma)\\
=&D_0+\sum_{h=1}^{+\infty}\,\sum_{j\in\Bscr_h } \frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\, ,
\end{aligned}
\end{equation}

where we define $j=h-2k$ and the set

\Bscr_h=\left\lbrace j\in\interi \, : \, j=
\begin{cases}
2n \quad \textit{if $h$ is even}\,, \, n\in\interi\,,\,-\frac{h}{2}\leq n \leq \frac{h}{2}\\
2n+1 \quad \textit{if $h$ is odd}\,, \, n\in\interi\,,\,\lfloor-\frac{h}{2}\rfloor\leq n \leq \lfloor\frac{h}{2}\rfloor
\end{cases} \right\rbrace\, , \qquad h\geq 1\, .

On the other hand, the previous set of indexes can be thought also as

\Bscr_h=
\begin{cases}
\Bscr_{h_-} \cup \Bscr_{h_+} \cup \lbrace j=0 \rbrace \quad \textit{if $h$ is even}\\
\Bscr_{h_-} \cup \Bscr_{h_+} \quad \textit{if $h$ is odd}
\end{cases} \, ,

where $\Bscr_{h_-}=\Bscr_{h}\cap \interi_{<0}$ and $\Bscr_{h_+}=\Bscr_{h}\cap \interi_{>0}\, $.
With the following notation, the sum appearing in Eq. (<ref>) can be decomposed as follows:

\begin{align*}
&\sum_{h=1}^{+\infty}\,\sum_{j\in\Bscr_h }=\sum_{\substack{h=1\\ h \textit{odd}}}^{+\infty}\,\left(\sum_{j\in \Bscr_{h_-}} + \sum_{j\in \Bscr_{h_+}}\right)+\sum_{\substack{h=2\\ h \textit{even}}}^{+\infty}\,\left(\sum_{j\in \Bscr_{h_-}} + \sum_{j\in \Bscr_{h_+}}+\sum_{j=0}\right)\, .
\end{align*}

Thus, observing that

\sum_{j\in \Bscr_{h_-}} \frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)=\sum_{j\in \Bscr_{h_+}} \frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\, ,

we can write Eq. (<ref>) as:

\begin{align*}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}=&\sum_{\substack{h=0\\h \textit{even}}}^{+\infty}\,\frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h}{2}\right)!\left(\frac{h}{2}\right)!}+\sum_{\substack{h=2\\h \textit{even}}}^{+\infty}\,\sum_{j\in\Bscr_{h_+} } \frac{D_{h}}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\notag\\
&+\sum_{\substack{h=1\\h \textit{odd}}}^{+\infty}\,\sum_{j\in\Bscr_{h_+} }\frac{D_h}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\, .
\end{align*}

Finally, reversing the order of the sums (over $j$ and $h\,$), we arrive at

\begin{equation}
\label{final1}
\begin{aligned}
\frac{1}{\left(1+\alpha^2-2\alpha \cos(\sigma)\right)^{\frac{l}{2}}}=&\sum_{\substack{h=0\\h \textit{even}}}^{+\infty}\,\frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h}{2}\right)!\left(\frac{h}{2}\right)!}+\sum_{\substack{j=2\\j \textit{even}}}^{+\infty}\,\sum_{h\in\Ascr_{j} } \frac{D_{h}}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\\
&+\sum_{\substack{j=1\\j \textit{odd}}}^{+\infty}\,\sum_{h\in\Ascr_{j} }\frac{D_h}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\, , \\
=&\sum_{\substack{h=0\\h \textit{even}}}^{+\infty}\,\frac{D_{h}}{2^h}\frac{h!}{\left(\frac{h}{2}\right)!\left(\frac{h}{2}\right)!}+\sum_{j=1}^{+\infty}\,\sum_{h\in\Ascr_{j} } \frac{D_{h}}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\cos(j\sigma)\, ,
\end{aligned}
\end{equation}

\Ascr_j=\left\lbrace h\in\naturali \, : \, h=
\begin{cases}
2i \quad \textit{if $j$ is even}\,, \, i\in\naturali\,,\,i\geq  \frac{j}{2}\\
2i+1 \quad \textit{if $j$ is odd}\,, \, i\in\naturali\,,\,i\geq \lfloor \frac{j}{2}\rfloor
\end{cases} \right\rbrace\, , \qquad j\geq 1\, .

Now we can finally put Eq. (<ref>) in the expression (<ref>) and compare the obtained quantity with (<ref>), having

\begin{align}
\label{coeff.Lapl.new}
&b_{s+\frac{1}{2}}^{(0)}(\alpha)=D_0+\sum_{h\in\Ascr_2}\frac{D_h}{2^{h}}\frac{h!}{\left(\frac{h}{2}\right)!\left(\frac{h}{2}\right)!}\, , & &b_{s+\frac{1}{2}}^{(j)}(\alpha)=\sum_{h\in\Ascr_j}\frac{D_h}{2^{h-1}}\frac{h!}{\left(\frac{h-j}{2}\right)!\left(\frac{h+j}{2}\right)!}\, ,
\end{align}

with $\Ascr_2=\left\lbrace h\in\naturali\, : \, h=2 i\,, \, i\in\naturali\,,\,i\geq 1 \right\rbrace\, $.
Finally, substituting in the previous expression (<ref>) the definition of $D_h$ (given by (<ref>)) and remembering that $l=2s+1\,$ and $\sigma=\lambda_2-\lambda_3\,$, we obtain (<ref>). This concludes the proof.

§ NORMAL FORM AROUND THE PERIODIC ORBITS A OR B

The normal form function $Z_{sec}^{(r)}$ and the normalizing transformation $\Phi^{(r)}$ of Eqs. (<ref>) and (<ref>) are computed by the following steps:

* Development of the Hamiltonian around the apsidal solution of the integrable part $\Hscr_{int}$: Let $(\psi_{\ast},\,\Gamma_{\ast},\,J_{\ast})$ be one of the two fixed points (apsidal corotation solutions) of the integrable Hamiltonian $\Hscr_{int}(\psi,\,\Gamma\,;J)\,$. We consider the translations

\begin{equation}
\label{trasla.apsidal.integ}
\psi=\psi_{\ast}+\varepsilon_*\delta\psi\, , \qquad \Gamma=\Gamma_{\ast}+\varepsilon_*\,\delta\Gamma\,,\qquad
J=J_{\ast}+\varepsilon_{J_*}\,\delta J\, ,
\end{equation}

where $\varepsilon_*$, $\varepsilon_{J_*}$ are both book-keeping symbols with numerical value equal to one. We then perform the following algebraic operations:

* substitute (<ref>) into the Hamiltonian $\Hscr_{sec}(\psi,\varphi,\Gamma,J)$ and expand the Hamiltonian in powers of the symbols $\varepsilon_*$, $\varepsilon_{J_*}$;

* in the resulting expression, replace $\cos(\varphi)$ and $\sin(\varphi)$ with $\varepsilon_*\cos(\varphi)$ and $\varepsilon_*\sin(\varphi)$ respectively;

* assign a unique book-keeping symbol to the expanded Hamiltonian $\Hscr_{sec}(\delta\psi,\,\varphi,\,\delta\Gamma,\,\delta J\,)$, by setting $\varepsilon_{J_*}^{a}\rightsquigarrow\,\lambda^{a-1}\,\, a\geq 1\,$, $\varepsilon_*\rightsquigarrow\,\lambda\,$;

* truncate the resulting expression at a maximum order $N_t$ in the book-keeping symbol $\lambda$.

After performing the above algebraic steps, the Hamiltonian resumes the form

\begin{equation}\label{ham0}
\begin{aligned}
\Hscr^{(0)}=\Hscr_{sec} &= c_1\,+Z_0
+\sum_{s=1}^{N_t}\lambda^s h_{s}^{(0)}(\delta\psi,\,\varphi,\,\delta\Gamma,\,\delta J\,) \\
&= c_1\,+\nu_*^{(0)}\delta J +\sum_{s=1}^{N_t}\lambda^s
\left(\sum_{m,\,k,\,n,\,l}\theta^{(0)}_{m,k,n,l}\,\delta\psi^{m}\delta\Gamma^{n} \delta J^{l}\,e^{i\,k\,\varphi} \right)
\end{aligned}
\end{equation}

where $c_1$ is a constant term, $\theta^{(0)}_{m,k,n,l}$ are constant coefficients, and $\nu_*^{(0)}<0$ is the (also constant) unperturbed frequency of the apsidal periodic orbit in the integrable approximation.

* Normal form to eliminate $\varphi$: The Hamiltonian (<ref>) is normalized by iterative steps using the method of composition of Lie series. The normalization algorithm is defined recursively, for $r=1,2,...$ by the relations:

\begin{equation}\label{hamrm1}
\Hscr^{(r-1)}=c_1\,+Z_0
+\sum_{s=1}^{r-1}\lambda^s Z_s(\delta\psi,\delta\Gamma,\delta J)
\left(\sum_{m,\,k,\,n,\,l}
\theta_{m,k,n,l}^{(r-1)}\,\delta\psi^{m}\delta\Gamma^{n} \delta J^{l}\,e^{i\,k\,\varphi} \right) ,
\end{equation}

\begin{equation}\label{chir}
\chi^{(r)}=\lambda^{r}\!\!\!\sum_{\substack{m,\,k,\,n,\,l\\ k\neq 0}}\frac{\theta^{(r-1)}_{m,k,n,l}}{i\,\nu_*^{(0)}\,k}\,
\delta\psi^{m}\delta\Gamma^{n} \delta J^{l}\,e^{i\,k\,\varphi}\, ,
\end{equation}

\begin{equation}\label{hamr}
\Hscr^{(r)}=\Bigg[\exp(L_{\chi^{(r)}})\Hscr^{(r-1)}\Bigg]^{\leq N_t} \, ,
\end{equation}

where $L_{\chi^{(r)}}$ denotes the Poisson bracket operator $L_{\chi^{(r)}}\cdot=\{\cdot,\chi^{(r)}\}$, and $[\cdot]^{\leq N_t}$ means truncation at the order $N_t$ in the book-keeping parameter $\lambda$.

* Final Hamiltonian and normalizing transformation: the normalizing transformation $\Phi^{(r)}$ is defined by:

\begin{eqnarray}\label{phifinal}
(\delta\psi,\varphi,\delta\Gamma,\delta J)&=&
\Phi^{(r)}(\delta\t\psi,\t\varphi,\delta\t\Gamma,\delta\t J) \\
\Bigg[
\left(\exp\left(L_{\chi^{(r)}}\right)\exp\left(L_{\chi^{(r-1)}}\right)\ldots\exp\left(L_{\chi^{(1)}}\right)(\delta\psi,\varphi,\delta\Gamma,\delta J)\right)
\Bigg|_{\substack{
\delta\psi={\delta\t\psi}\\
\varphi=\t{\varphi}\\
\delta\Gamma={\delta\t\Gamma}\\
\delta J={\delta\t{J}}}}
\Bigg]^{\leq N_t}~ \nonumber
\end{eqnarray}

and the final Hamiltonian

\begin{equation}\label{hamfinal}
\begin{aligned}
\Hscr^{(r)}&=
\Bigg[
\left(\exp\left(L_{\chi^{(r)}}\right)\exp\left(L_{\chi^{(r-1)}}\right)\ldots\exp\left(L_{\chi^{(1)}}\right)
\Hscr^{(0)}\right)
\Bigg|_{\substack{
\delta\psi={\delta\t\psi}\\
\varphi=\t{\varphi}\\
\delta\Gamma={\delta\t\Gamma}\\
\delta J={\delta\t{J}}}}
\Bigg]^{\leq N_t}~ \\
+R^{(r)}(\delta\t{\psi},\,\t{\varphi},\,\delta\t{\Gamma},\,\delta\t{J}\,)\, ,
\end{aligned}
\end{equation}

where $1\leq r\leq N_t$ is chosen so that the normal form term $Z^{(r)}(\delta\t{\psi},\,\delta\t{\Gamma},\,\delta\t{J}\,)$
yields the best possible approximation to the dynamics.

Note that the constant term $c_1$ produced at step 1 above has to be carried along all the successive Hamiltonians $\Hscr^{(r)}$ produced through the normalization process, being eventually included in the normal form $Z^{(r)}$. This is necessary, since this constant appears in the third of the algebraic equations (<ref>) by which the fixed point corresponding to the apsidal periodic orbit is computed in the new canonical variables.

§ PROOF OF PROPOSITION <REF>

The eigenvalues of the matrix M (<ref>) are $\lambda_{1,2}=\pm \sqrt{\mathfrak{a}\mathfrak{b}}\,$; then if $\mathfrak{a}\mathfrak{b}>0\,$, then the eigenvalues are real and opposite, instead if $\mathfrak{a}\mathfrak{b}<0\,$, then the eigenvalues are complex and conjugate. We have to analyze all the possible cases. We can have

\begin{align}
\mathfrak{a}>0\\
\mathfrak{b}>0
\end{cases}\, , &
\mathfrak{a}<0\\
\mathfrak{b}<0
\end{cases} \quad \Longrightarrow\quad \lambda_{1,2}\quad \mathit{reals} \label{autoval.real}
\\
\mathfrak{a}>0\\
\mathfrak{b}<0
\end{cases}\, , &
\mathfrak{a}<0\\
\mathfrak{b}>0
\end{cases} \quad \Longrightarrow\quad \lambda_{1,2}\quad \mathit{immaginary}\, . \label{autoval.imm}
\end{align}

Let us start to understand the change of sign of $\mathfrak{b}\,$; remembering the definition (<ref>)

\begin{equation*}
\mathfrak{b}=\frac{3 \mathcal{G}^2 L_2^3\, (3 G_3^2 - L_2^2 + L_z^2)\, m_0\, m_3^7}{8 \,G_3^5\, L_3^3\, m_2^3}>0 \quad\Longleftrightarrow\quad  \tilde{\mathfrak{b}}:= L_z^2 +3L_3^2(1-\e_3^2)-L_2^2>0\, ,
\end{equation*}

being $(3 \mathcal{G}^2 L_2^3 m_0 m_3^7)/(8 (1-\e_3^2)^{5/2} L_3^8 m_2^3)>0\,$. Now, let us observe that if

\begin{align}
\label{segnob1}
C:=L_2^2-3L_3^2(1-\e_3^2)<0\quad\Longrightarrow\quad \e_3^2<1-\frac{L_2^2}{3 L_3^2}\, ,
\end{align}

then $\tilde{\mathfrak{b}}>0$ automatically; in order to have that, it is necessary that $1-\frac{L_2^2}{3 L_3^2}>0\,$, i.e. to have $0< L_2< \sqrt{3}L_3\,$. Thus, we can conclude that $\mathfrak{b}>0$ ( or equivalently $\tilde{\mathfrak{b}}>0\,$) if one of the two is fullfilled

i) $0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,0\leq \e_3\leq\sqrt{1-\frac{L_2^2}{3L_3^2}}\,$, with $L_z^2>0\,$ or

ii) $\left(0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,\e_3>\sqrt{1-\frac{L_2^2}{3L_3^2}}\right)\,\,\vee\,\,L_2\geq\sqrt{3}L_3\,$, with $L_z^2>C\,$.

Now, let us study the sign of $\mathfrak{a}\,$; as before, from definition (<ref>)

\begin{equation}
\begin{split}
&\mathfrak{a}=-\frac{3 \mathcal{G}^2 L_2 (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) m_0 m_3^7}{16 G_3^5 L_3^3 m_2^3}>0\quad\Longleftrightarrow\\
&\tilde{\mathfrak{a}}:= L_3^4(1-\e_3^2)^2 - \frac{4}{5} L_3^2 L_2^2(1-\e_3^2) + \frac{3}{5} L_2^4 - 2 L_3^2(1-\e_3^2) L_z^2 -
\frac{8}{5} L_2^2 L_z^2 +  L_z^4<0\, ,
\end{split}
\end{equation}

being $-(3 \mathcal{G}^2 L_2\, m_0 \,m_3^7)/(16\, L_3^8(1-\e_3)^{5/2} m_2^3)<0\,$. Then $\mathfrak{a}>0$ (or equivalently $\tilde{\mathfrak{a}}<0\,$) iff $A<L_z^2<B\,$,
where $A$ and $B$ are defined as in (<ref>), i.e.

\begin{align*}
&A=\frac{1}{5}\left(4L_2^2+5 L_3^2(1-\e_3^2)-L_2\sqrt{L_2^2+60 L_3^2(1-\e_3^2)}\right)\, ,\\
&B=\frac{1}{5}\left(4L_2^2+5 L_3^2(1-\e_3^2)+L_2\sqrt{L_2^2+60 L_3^2(1-\e_3^2)}\right)\, .
\end{align*}
Let us observe that $A$ and $B$ are both greater than zero. Infact, being $0\leq \e_3 < 1\,$, it is obvious that $B>0\,$; on the other hand $A>0$ iff

\begin{align*}
&\sqrt{L_2^2+60 L_3^2(1-\e_3^2)}<4L_2+\frac{5 L_3^2}{L_2}(1-\e_3^2) & &\Longleftrightarrow\\
&L_2^2+60 L_3^2(1-\e_3^2)<16 L_2^2 + 25 \frac{L_3^4}{L_2^2}(1-\e_3^2)^2+40 L_3^2(1-\e_3^2) & &\Longleftrightarrow \\
& 3L_2^4-4 L_3^2 (1-\e_3^2) L_2^2 + 5 (1-\e_3^2)^2 L_3^4 > 0
\end{align*}

that holds $\forall\,\, L_2\, ,L_3\,, 0\leq \e_3< 1\,$, being $\Delta=-44 L_3^4 (1-\e_3^2)^2 <0\,$. Then, we can conclude that if we are in the case i), (i.e. $\displaystyle{0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,0\leq \e_3\leq\sqrt{1-\frac{L_2^2}{3L_3^2}}}\,$) and $A< L_z^2< B\,$, then $\mathfrak{b}>0\,$ and $\mathfrak{a}>0\,$, then we have reals eigenvalues; instead if, in the same case i), $0<L_z^2< A\,$ or $L_z^2> B\,$, then $\mathfrak{b}>0\,$, $\mathfrak{a}<0\,$, having complex immaginary eigenvalues. This proves the Case i) of the Proposition.

In the case ii), i.e. $\left(0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,\e_3>\sqrt{1-\frac{L_2^2}{3L_3^2}}\right)\,\,\vee\,\,L_2\geq\sqrt{3}L_3\,$, we have to solve (<ref>) and (<ref>), that become, respectively

\begin{align}
\mathfrak{a}>0\\
\mathfrak{b}>0
\end{cases}
\, =\,
\begin{cases}
A<L_z^2< B\\
\end{cases}\, , &
\mathfrak{a}<0\\
\mathfrak{b}<0
\end{cases} \, =\,
\begin{cases}
0<L_z^2< A \, \vee\, L_z^2>B\\
\end{cases}\, ,\label{autoval.real2}
\\
\mathfrak{a}>0\\
\mathfrak{b}<0
\end{cases}\, =\,
\begin{cases}
A<L_z^2< B\\
\end{cases}\, , &
\mathfrak{a}<0\\
\mathfrak{b}>0
\end{cases}\, =\,
\begin{cases}
0<L_z^2< A \, \vee\, L_z^2>B\\
\end{cases}\, ;\label{autoval.imm2}
\end{align}

then, we have to understand the position of the value $C$ with respect to $A$ and $B\,$.

Case 1 : $C< A$

$$C< A\quad\Longleftrightarrow\quad \sqrt{L_2^2+60 L_3^2(1-\e_3^2)}< -L_2+ \frac{20 L_3^2}{L_2}(1-\e_3^2)$$

\begin{align*}
\begin{cases}
\displaystyle{-L_2+ \frac{20 L_3^2}{L_2}(1-\e_3^2)\geq 0}\\
\displaystyle{L_2^2+60 L_3^2(1-\e_3^2)< L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,
\end{cases}
\quad\Longleftrightarrow\quad
\begin{cases}
\displaystyle{\e_3^2\leq 1-\frac{L_2^2}{20 L_3^2} }\\
\displaystyle{\e_3^2< 1-\frac{L_2^2}{4 L_3^2} }\,,
\end{cases}
\end{align*}

then if $\e_3^2 < 1-\frac{L_2^2}{4 L_3^2}\, $. The last inequality makes sense if $1-\frac{L_2^2}{4 L_3^2}> 0\,$, that is $0\leq L_2< 2 L_3\,$. Remembering that we are in the case ii), it follows that

\begin{align*}
C<A \quad\Longleftrightarrow\quad &\left(0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,\sqrt{1-\frac{L_2^2}{3L_3^2}}< \e_3 < \sqrt{1-\frac{L_2^2}{ 4 L_3^2}} \right)\, \, \vee \, \,\\
& \left(\sqrt{3}L_3\leq L_2 < 2 L_3 \,\,\wedge\,\,0\leq \e_3 < \sqrt{1-\frac{L_2^2}{ 4 L_3^2}}\right)\, .
\end{align*}

Therefore, studying (<ref>) and (<ref>), we can conclude that if
$\left(0< L_2<\sqrt{3}L_3 \,\,\wedge\,\,\sqrt{1-\frac{L_2^2}{3L_3^2}}< \e_3 < \sqrt{1-\frac{L_2^2}{ 4 L_3^2}} \right)$  $\vee$ $ \left(\sqrt{3}L_3\leq L_2 < 2 L_3 \,\,\wedge\,\,0\leq \e_3 < \sqrt{1-\frac{L_2^2}{ 4 L_3^2}}\right)\,$,
then the eigenvalues are reals and opposite if $0<L_z^2<C\,$ (see the second of (<ref>)), purely immaginary and coniugates if $C<L_z^2<A\,$ (see the second of (<ref>)), reals and opposite if $A<L_z^2<B$ (see the first of (<ref>)) and again purely immaginary and coniugates if $L_z^2>B\,$ (see the second of (<ref>)), proving the Case ii) of the Proposition.

Case 2 : $C= A$

Following the same calculations done for the previous case ($C<A\,$), it easily follows that

\begin{align*}
C= A\quad\Longleftrightarrow\quad &\left(0< L_2\leq 2 L_3 \,\,\wedge\,\, \e_3 = \sqrt{1-\frac{L_2^2}{ 4 L_3^2}} \right) \,\, ;
\end{align*}

in this case, studying (<ref>) and (<ref>), we can conclude that the eigenvalues are reals and opposite if $0<L_z^2<C=A\,$ (see the second of (<ref>)), again reals and opposite if $C=A<L_z^2<B$ (see the first of (<ref>)) (i.e. are reals and opposite if $0<L_z^2<B\,$, passing from $0$ if $L_z^2=A=C\,$) and purely immaginary and coniugates if $L_z^2>B\,$ (see the second of (<ref>)). This proves the Case iii) of the Proposition.

Case 3 : $A< C< B$

\begin{align}
\label{prima}
&C> A\quad&\Longleftrightarrow&\quad \sqrt{L_2^2+60 L_3^2(1-\e_3^2)}> -L_2+ \frac{20 L_3^2}{L_2}(1-\e_3^2)\\
\label{seconda}
&C< B\quad&\Longleftrightarrow&\quad \sqrt{L_2^2+60 L_3^2(1-\e_3^2)}> L_2- \frac{20 L_3^2}{L_2}(1-\e_3^2)\, .
\end{align}

Let us start from (<ref>); we have to solve the following

\begin{align*}
\displaystyle{-L_2+ \frac{20 L_3^2}{L_2}(1-\e_3^2)\geq 0}\\
\displaystyle{L_2^2+60 L_3^2(1-\e_3^2)> L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,
\end{cases}
\quad \wedge\quad
\begin{cases}
\displaystyle{-L_2+ \frac{20 L_3^2}{L_2}(1-\e_3^2)< 0}\\
\phantom{\displaystyle{L_2^2+60 L_3^2(1-e_3^2)> L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,}
\end{cases}
\end{align*}

that gives

\begin{equation}
\label{e3.prima}
\left(\begin{cases}
\displaystyle{\e_3^2\leq 1-\frac{L_2^2}{20 L_3^2} }\\
\displaystyle{\e_3^2> 1-\frac{L_2^2}{4 L_3^2} }\,
\end{cases}
\quad \wedge\quad
\begin{cases}
\displaystyle{\e_3^2> 1-\frac{L_2^2}{20 L_3^2} }\\
\phantom{\displaystyle{\e_3^2> 1-\frac{L_2^2}{4 L_3^2} }\,,}
\end{cases}\right)\qquad\Longrightarrow\qquad \e_3^2> 1-\frac{L_2^2}{4 L_3^2}\, .
\end{equation}

Similarly, we can solve (<ref>), having the following

\begin{align*}
\displaystyle{L_2- \frac{20 L_3^2}{L_2}(1-\e_3^2)\geq 0}\\
\displaystyle{L_2^2+60 L_3^2(1-\e_3^2)> L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,
\end{cases}
\quad \wedge\quad
\begin{cases}
\displaystyle{L_2- \frac{20 L_3^2}{L_2}(1-\e_3^2)< 0}\\
\phantom{\displaystyle{L_2^2+60 L_3^2(1-\e_3^2)> L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,}
\end{cases}
\end{align*}

that gives

\begin{equation}
\label{e3.seconda}
\left(\begin{cases}
\displaystyle{\e_3^2\geq 1-\frac{L_2^2}{20 L_3^2} }\\
\displaystyle{\e_3^2> 1-\frac{L_2^2}{4 L_3^2} }\,
\end{cases}
\quad \wedge\quad
\begin{cases}
\displaystyle{\e_3^2< 1-\frac{L_2^2}{20 L_3^2} }\\
\phantom{\displaystyle{\e_3^2> 1-\frac{L_2^2}{4 L_3^2} }\,,}
\end{cases}\right)\qquad\Longrightarrow\qquad \forall\,\e_3 .
\end{equation}

Finally, putting together (<ref>) and (<ref>), we get $\e_3^2> 1-\frac{L_2^2}{4 L_3^2}\,$. Let us observe that if $1-\frac{L_2^2}{4 L_3^2}< 0\,$, i.e. $ L_2> 2 L_3\,$,
then the last inequality for $\e_3$ is automatically satisfy. Then, remembering that we are in the case ii) and putting togheter the results, it follows that

\begin{align*}
A<C<B \quad\Longleftrightarrow\quad &\left(0< L_2\leq 2 L_3 \,\,\wedge\,\, \e_3 > \sqrt{1-\frac{L_2^2}{ 4 L_3^2}} \right)\, \, \vee \, \,\\
& \quad L_2 > 2 L_3\, .
\end{align*}

Again, studying (<ref>) and (<ref>), we can conclude that if $\displaystyle{\left(0< L_2\leq 2 L_3 \,\,\wedge\,\, \e_3 > \sqrt{1-\frac{L_2^2}{ 4 L_3^2}} \right)}$  $\vee$ $ L_2 > 2 L_3\,$, then the eigenvalues are reals and opposite if $0<L_z^2<A\,$ (see the second of (<ref>)), purely immaginary and coniugates if $A<L_z^2<C\,$ (see the first of (<ref>)), reals and opposite if $C<L_z^2<B$ (see the first of (<ref>)) and again purely immaginary and coniugates if $L_z^2>B\,$ (see the second of (<ref>)). This proves the last Case iv) of the Proposition.

Actually, there is a last case to be analyzed, that will not give any contribution, i.e.

Case 4 : $ C\geq  B$

$$C\geq B\quad\Longleftrightarrow\quad \sqrt{L_2^2+60 L_3^2(1-\e_3^2)}\leq L_2- \frac{20 L_3^2}{L_2}(1-\e_3^2)$$

\begin{align*}
\begin{cases}
\displaystyle{L_2- \frac{20 L_3^2}{L_2}(1-\e_3^2)\geq 0}\\
\displaystyle{L_2^2+60 L_3^2(1-\e_3^2)\leq L_2^2+ \frac{400 L_3^4}{L_2^2}(1-\e_3^2)^2- 40 L_3^2 (1-\e_3^2)}\,
\end{cases}
\quad\Longleftrightarrow\quad
\begin{cases}
\displaystyle{\e_3^2\geq 1-\frac{L_2^2}{20 L_3^2} }\\
\displaystyle{\e_3^2\leq 1-\frac{L_2^2}{4 L_3^2} }
\end{cases}
\end{align*}

that has no solutions. This concludes the proof of the Proposition.