---
abstract: 'The Lidov-Kozai mechanism allows a body to periodically exchange its eccentricity with inclination. It was first discussed in the framework of the quadrupolar secular restricted three-body problem, where the massless particle is the inner body, and later extended to the quadrupolar secular nonrestricted three body problem [@1969CeMec...1..200H; @1976CeMec..13..471L; @1994CeMDA..58..245F]. In this paper, we propose a different point of view on the problem by looking first at the restricted problem where the massless particle is the outer body. In this situation, equilibria at high mutual inclination appear [[@2006PhyD..213...15P]]{}, which correspond to the population of stable particles that @2008MNRAS.390.1377V [@2009MNRAS.394.1721V] find in stable, high inclination circumbinary orbits around one of the components of the quadruple star HD 98800. We provide [a simple analytical framework using a vectorial formalism for these situations. We also look]{} at the evolution of these high inclination equilibria in the non restricted case.'
author:
- |
    F. Farago[^1] and J. Laskar\
    ASD, IMCCE-CNRS UMR8028, Observatoire de Paris, UPMC, 77 avenue Denfert-Rochereau, 75014 Paris, France
bibliography:
- 'FL2009.bib'
title: 'High inclination orbits in the secular quadrupolar three-body problem'
---

celestial mechanics – planetary systems – methods: analytical – methods: $N$-body simulations.

Introduction {#se:intro}
============

As it is known, the secular three-body problem after node reduction has two degrees of freedom (e. g. [@1905QB351.P73......; @2002CeMDA..84..283M]). However, due to what @1976CeMec..13..471L called a [*happy coincidence*]{}, this problem is integrable when it is expanded up to order two in the ratio of semi-major axes, i.e. at the quadrupolar approximation. Indeed, the argument of perihelion of the outer body does not explicitly appear in the quadrupolar expansion of the secular problem, thus giving one more integral of motion linked to the eccentricity of the outer body.

The limiting case where the inner body has no mass has been extensively studied . We will call this problem the [*inner restricted problem*]{}, while the converse case where the two inner bodies are massive and the outer body is massless will be called the [*outer restricted problem*]{}. In the inner restricted case, the conservation of the normal component of the angular momentum enables the inner particle to periodically exchange its eccentricity with inclination (the so-called Lidov-Kozai mechanism). The inner restricted model is well suited when the inner body has a small mass with respect to the other two. However, when looking at higher mass ratios, for example in triple star systems, this is no longer justified.

Since the Hamiltonian of the quadrupolar problem of three masses is very similar to that of the inner restricted problem when it is written in elliptic variables, the study of the massive problem has mainly focused on the dynamics of the two inner bodies [@1969CeMec...1..200H; @1976CeMec..13..471L; @1994CeMDA..58..245F]. These previous works completely classified the different dynamical regimes and bifurcations, using the equations of motion of the inner binary.

There is however another limit-case to the massive problem, which is the outer restricted problem. [@2006PhyD..213...15P have studied this case and discussed the existence and stability of equilibria in the non-averaged system using the framework of KAM theory. We give here a very simple model of the outer restricted case which provides an alternate formulation of these previous results and which is directly usable in an astronomical context. We also fully describe the possible motions of the bodies and give an analytical expression of their frequencies.]{} We use this model to explain the results of @2008MNRAS.390.1377V [@2009MNRAS.394.1721V], who find populations of particles at very high inclinations around one of the components of the double-binary star HD 98800, which are stable even under the perturbation of the other component. We then look at the quadrupolar problem of three masses from the perspective of the outer restricted problem [and show how the inner and outer restricted cases are related to the general case]{}. Vectorial methods as developed by @2006Icar..185..312B [@2009Icar..201..750B; @2009AJ....137.3706T] are extremely well suited for this approach.

Secular outer restricted problem {#se:rtbp}
================================

Derivation of the Hamiltonian
-----------------------------

We consider here the case of a massless particle orbiting a central binary object. We do not restrict ourselves with respect to inclinations or eccentricities. The components of the binary have masses $m_0$ and $m_1$, the binary’s total mass is ${\mathcal{M}_{01}}=m_0+m_1$ and its reduced mass is ${\beta_1}=m_0m_1/(m_0+m_1)$. The two massive components have barycentric positions ${\mathbf{u}}_0\ \text{and}\ {\mathbf{u}}_1$. We also denote $\delta=m_0/{\mathcal{M}_{01}}$ and ${{\mathbf{r}}_1}={\mathbf{u}}_1-{\mathbf{u}}_0$, and ${\mathbf{r}}_2$ is the position of the outer particle relatively to the barycentre of the inner binary. Using these notations, the particle has the following Hamiltonian:

$$H = \frac{\tilde{{\mathbf{r}}}^2_2}{2}-G\left(\frac{m_0}{|{\mathbf{r}}_2-{\mathbf{u}}_0|}+\frac{m_1}{|{\mathbf{r}}_2-{\mathbf{u}}_1|}\right)\ ,$$

where $\tilde{{\mathbf{r}}}_2=\dot{{\mathbf{r}}}_2$ is the canonical momentum associated to the position of the massless particle, ${\mathbf{r}}_2$. Since ${\mathbf{u}}_0=-(1-\delta){{\mathbf{r}}_1}$ and ${\mathbf{u}}_1=\delta{{\mathbf{r}}_1}$, we can rewrite the Hamiltonian as:

$$H=\frac{\tilde{{\mathbf{r}}}^2_2}{2}-G\left(\frac{m_0}{|{\mathbf{r}}_2+(1-\delta){{\mathbf{r}}_1}|}
 +\frac{m_1}{|{\mathbf{r}}_2-\delta{{\mathbf{r}}_1}|}\right)\ .$$

We now suppose that ${r_1}\ll r_2$ and expand the Hamiltonian to order 2 in $r_1/r_2$:

$$\label{eq:hamqp}
 H = \left(\frac{\tilde{{\mathbf{r}}}^2_2}{2}-\frac{G{\mathcal{M}_{01}}}{r_2}\right)
 -\frac{G{\beta_1}}{2r_2^3}\left(3\frac{{\mathopen ( {\mathbf{r}}_2 . {{\mathbf{r}}_1}\mathclose )}^2}{r_2^2}-{r_1}^2\right)\ .$$

The first term is the Keplerian energy of the particle interacting with the binary, seen as a point mass ${\mathcal{M}_{01}}$. It is equal to $-G{\mathcal{M}_{01}}/2a_2$, where $a_2$ is the semi major axis of the particle.

Since we are interested in the secular behaviour of the particle, we average this quadrupolar Hamiltonian over the mean anomalies of the binary (${M_1}$) and of the particle ($M_2$). In order to do this, we first introduce four unit vectors: $({\mathbf{i}},{\mathbf{j}},{\mathbf{k}})$ are bound to the orbit of the binary, remain constant, and will provide a natural reference frame; ${\mathbf{w}}$ is bound to the orbit of the particle and will vary. More precisely, ${\mathbf{i}}$ points in the direction of the perihelion of the binary, ${\mathbf{k}}$ is colinear to the angular momentum of the binary, and ${\mathbf{j}}={{\mathbf{k}} \wedge {\mathbf{i}}}$; the last vector ${\mathbf{w}}$ is colinear to the angular momentum of the massless particle.

We can then compute the following averaged quantities, where quantities indexed with $1$ relate to the binary, quantities with index $2$ relate to the particle, and ${\mathbf{u}}$ is an arbitrary fixed vector (see for instance the appendix of @2006Icar..185..312B):

$$\begin{aligned}
 {\left\langle {r_1}^2 \right\rangle}_{{M_1}} & = & {a_1}^2\left(1+\frac{3}{2}{e_1}^2\right)\ ; \\
 {\left\langle {\mathopen ( {{\mathbf{r}}_1}. {\mathbf{r}}_2 \mathclose )}^2 \right\rangle}_{{M_1}} & = & \frac{{a_1}^2}{2}(r_2^2-{\mathopen ( {\mathbf{k}} . {\mathbf{r}}_2 \mathclose )}^2)\nonumber\\
    &&+\frac{{a_1}^2{e_1}^2}{2}(4{\mathopen ( {\mathbf{i}} . {\mathbf{r}}_2 \mathclose )}^2 - {\mathopen ( {\mathbf{j}} . {\mathbf{r}}_2 \mathclose )}^2)\ ; \\
 {\left\langle \frac{1}{r_2^3} \right\rangle}_{M_2} & = & \frac{1}{a_2^3(1-e_2^2)^{3/2}}\ ; \\
 {\left\langle \frac{{\mathopen ( {\mathbf{r}}_2 . {\mathbf{u}} \mathclose )}^2}{r_2^5} \right\rangle}_{M_2} & = & 
    \frac{u^2-{\mathopen ( {\mathbf{w}} . {\mathbf{u}} \mathclose )}^2}{2a_2^3(1-e_2^2)^{3/2}}\ .\end{aligned}$$

The substitution of these expressions in yields:

$$\begin{gathered}
\label{eq:hmoy}
{\left\langle H \right\rangle}_{{M_1},M_2} = -\frac{G{\mathcal{M}_{01}}}{2a_2} -\frac{3}{8}\frac{G{\beta_1}{a_1}^2}{a_2^3(1-e_2^2)^{3/2}}\times\\
\left[\left({e_1}^2-\frac{1}{3}\right)
+{\mathopen ( {\mathbf{k}} . {\mathbf{w}} \mathclose )}^2-{e_1}^2(4{\mathopen ( {\mathbf{i}} . {\mathbf{w}} \mathclose )}^2-{\mathopen ( {\mathbf{j}} . {\mathbf{w}} \mathclose )}^2)\right]\end{gathered}$$

Since the particle has no mass, the only variable element of the binary is its mean anomaly ${M_1}$. After averaging over this angle, it is no longer present in the Hamiltonian. After averaging over the mean anomaly of the particle, its semi major axis $a_2$ becomes constant. Moreover, ${\mathbf{w}}=\sin i_2 \sin\Omega_2\ {\mathbf{i}} -
\sin i_2 \cos\Omega_2\ {\mathbf{j}} + \cos i_2\ {\mathbf{k}},$ so the argument of pericentre $\omega_2$ of the particle does not appear in the averaged Hamiltonian. Hence, at the quadrupolar order, the conjugate momentum associated to $\omega_2$, i.e. the norm of the angular momentum of the particle $G_2=\sqrt{G{\mathcal{M}_{01}}a_2(1-e_2^2)}$, is constant. Therefore the eccentricity $e_2$ of the particle is constant. This fact is a feature of the quadrupolar expansion, not a property of the restricted problem. As such it remains true when the outer body has a non-zero mass (see section \[se:ftbp\]). This is the [*happy coincidence*]{} that @1976CeMec..13..471L noted. Finally, only one degree of freedom remains, related to the couple $(i_2,\Omega_2)$.

If we drop the constant terms in , and introduce the mean motion $n_1$ of the binary into the Hamiltonian ($n_1^2{a_1}^3=G{\mathcal{M}_{01}}$), we get the following expression[^2] [(see also eq. 10 in [@2006PhyD..213...15P])]{}:

$${\left\langle H \right\rangle} = -\frac{\alpha G_2}{2}
\left[{\mathopen ( {\mathbf{k}} . {\mathbf{w}} \mathclose )}^2-{e_1}^2(4{\mathopen ( {\mathbf{i}} . {\mathbf{w}} \mathclose )}^2-{\mathopen ( {\mathbf{j}} . {\mathbf{w}} \mathclose )}^2) \right]\ ,$$

where $$\alpha = \frac{3}{4}n_1 \left(\frac{{a_1}}{a_2}\right)^{7/2} \frac{{\beta_1}}{{\mathcal{M}_{01}}}\frac{1}{(1-e_2^2)^2}\ .$$

This Hamiltonian can be rewritten in a very compact form as: $${\left\langle H \right\rangle} = -\frac{1}{2}\phantom{.}^t{\mathbf{w}}.{\mathbf{T}}.{\mathbf{w}}\ ,\ \text{where:}$$ $${\mathbf{T}} = \alpha G_2\left(\begin{array}{ccc}
                -4{e_1}^2 & 0 & 0 \\
                0 & {e_1}^2 & 0 \\
                0 & 0 & 1
               \end{array}\right)\ .$$

We also give the expression of the Hamiltonian using the inclination and the node of the particle:

$$\label{eq:HrestIO}
 {\left\langle H \right\rangle}=-\frac{\alpha G_2}{4}\left[2\cos^2i_2 - {e_1}^2\sin^2i_2\left(3-5\cos2\Omega_2\right)\right]\ .$$

Equations of motion {#se:restreqns}
-------------------

As discussed in [@2006Icar..185..312B], the equations of motion for ${\mathbf{w}}$ are simply obtained by:

$$\dot{{\mathbf{w}}} = \frac{1}{G_2}\nabla_{{\mathbf{w}}}{\left\langle H \right\rangle}\wedge{\mathbf{w}}\ .$$

After computing the gradient, we find:

$$\dot{{\mathbf{w}}} = -\alpha
  \left[(\mathbf{k}.\mathbf{w})(\mathbf{k}\wedge\mathbf{w})
       -{e_1}^2(4(\mathbf{i}.\mathbf{w})(\mathbf{i}\wedge\mathbf{w})
              -(\mathbf{j}.\mathbf{w})(\mathbf{j}\wedge\mathbf{w}))\right]\ .$$

If we note $x=(\mathbf{i}.\mathbf{w})$, $y=(\mathbf{j}.\mathbf{w})$, and $z=(\mathbf{k}.\mathbf{w})$, we get the following system for $(x,y,z)$:

$$\begin{aligned}
\dot{x} & = & \alpha(1-{e_1}^2)yz\ ; \label{eq:xdot}\\
\dot{y} & = & -\alpha(1+4{e_1}^2)zx\ ; \\
\dot{z} & = & 5\alpha {e_1}^2 xy\ .\label{eq:zdot}\end{aligned}$$

In these variables, the fact that ${\mathbf{w}}$ is a unit vector and the energy integral translate into the following equalities:

$$\begin{aligned}
 x^2+y^2+z^2 & = & 1\ ; \label{eq:wsphere}\\
 z^2 - {e_1}^2(4x^2-y^2) & = & h = \text{Cst}\ . \label{eq:Hsurf}\end{aligned}$$

The system of three differential equations – has thus two independent first integrals and is as such integrable. It is also straightforward from these two relations that $$\label{eq:hineg}
    -4{e_1}^2\leqslant h\leqslant 1\ .$$

Motion of a massless body around a circular binary
--------------------------------------------------

In the case of a circular binary, the Hamiltonian and the equations of motion greatly simplify[^3]:

$${\left\langle H \right\rangle} = -\frac{\alpha G_2}{2}{\mathopen ( {\mathbf{k}} . {\mathbf{w}} \mathclose )}^2\ ,$$

$$\dot{{\mathbf{w}}} = -\alpha
  (\mathbf{k}.\mathbf{w})(\mathbf{k}\wedge\mathbf{w})$$

The scalar product $(\mathbf{k}.\mathbf{w})=\cos i_2$ remains constant, and the nodes of the orbit of the particle simply precess around the angular momentum of the binary, with a constant precession rate: $$\dot{\Omega}=-\alpha\cos i_2 =
 -\frac{3}{4}n_1 \left(\frac{{a_1}}{a_2}\right)^{7/2} \frac{{\beta_1}}{{\mathcal{M}_{01}}}\frac{\cos i_2}{(1-e_2^2)^2}$$

This precession is equivalent to the precession generated by the quadrupolar potential of a circular and homogeneous ring of mass $\beta_1$ and of radius $a_1$ [following an idea which can be traced back to Gauss (see [@2009MNRAS.394.1085T] and references therein)]{}.

Motion of a massless body around an elliptic binary
---------------------------------------------------

### Qualitative overview

When the binary is elliptic, the situation changes and cannot be explained any longer by the quadrupolar torque of a circular ring. If we substitute $z^2$ in using , we get: $$\begin{aligned}
\label{eq:hinterw}
 (1+4{e_1}^2)x^2+(1-{e_1}^2)y^2 &=1-h\geqslant 0\ ,\\
  x^2+y^2+z^2 &= 1\ .\end{aligned}$$

The intersections of the energy surfaces and the normalized angular momentum sphere of the particle can thus be seen as the intersections of elliptic cylinders with the unit sphere. For a given value of the energy $h$, the extremity of the unit angular momentum vector of the particle ${\mathbf{w}}$ will move on the intersection of the corresponding cylinder with the unit sphere. Figures \[fig:Hxy\].a and c show these intersections for different values of the energy as dotted lines drawn on the unit sphere, in two situations where the binary has an eccentricity of $0.5$ and $0.2$ respectively. The three axes correspond to the scalar products $x$, $y$ and $z$ that are defined in section \[se:restreqns\].

![\[fig:Hxy\]Intersections of the energy surfaces and the unit angular momentum sphere (a) and its projection in the $(x,y)$ plane (b) for ${e_1}=0.5$. Intersections of the energy surfaces and the unit angular momentum sphere (c) and its projection in the $(x,y)$ plane (d) for ${e_1}=0.2$.](FL2009-fig1){width="\columnwidth"}

There are four visible kinds of trajectories: closed trajectories around the two poles of the sphere $(x,y,z)=(0,0,\pm 1)$, and closed trajectories around the points $(x,y,z)=(\pm 1,0,0)$.

When the extremity of the angular momentum of the particle ${\mathbf{w}}$ follows a trajectory around the north pole, it means that it precesses around the angular momentum of the binary ${\mathbf{k}}$ with an inclination that is strictly inferior to $90^\circ$: in this case, the orbital motion of the particle is prograde relatively to the orbital motion of the binary.

When the extremity of the angular momentum of the particle ${\mathbf{w}}$ follows a trajectory around the south pole, it means that it precesses around the opposite of the angular momentum of the binary, $-{\mathbf{k}}$, with an inclination that is strictly superior to $90^\circ$: in this case, the orbital motion of the particle is retrograde relatively to the orbital motion of the binary.

When the extremity of the angular momentum of the particle ${\mathbf{w}}$ follows a trajectory around one of the two points $(x,y,z)=(\pm 1,0,0)$, it precesses around the direction of the perihelion of the binary or the opposite of this direction. In this case, the inclination oscillates around $\pm 90^\circ$.

### Frequencies

The frequencies of these motions can be found analytically. Indeed, using equation , we see that $x$ and $y$ are on ellipses or arcs of ellipses bounded by the unit circle (figures \[fig:Hxy\].b and d show respectively the cases where ${e_1}=0.5$ and $0.2$). Thus, there is an angle $\phi$ such that:

$$\begin{aligned}
 x & = & \sqrt{\frac{1-h}{1+4{e_1}^2}} \cos\phi\ , \label{eq:xell}\\
 y & = & \sqrt{\frac{1-h}{1-{e_1}^2}} \sin \phi\ .\end{aligned}$$

Using , we can then express $z^2$ as:

$$\label{eq:z2phi}
 z^2 = \frac{h+4{e_1}^2}{1+4{e_1}^2} - \frac{5{e_1}^2}{1+4{e_1}^2}\frac{1-h}{1-{e_1}^2} \sin^2\phi\ .$$

There are two opposite values of $z$ for each $\phi$, reflecting the symmetry of the system with respect to the orbital plane of the binary. If we use expression in combination with equation , we obtain the following equation for $\dot{\phi}$:

$$\begin{aligned}
 \dot{\phi}&=&\mp\alpha\sqrt{(1-{e_1}^2)(h+4{e_1}^2)}
\sqrt{1-k^2\sin^2\phi}\ , \label{eq:dotphi}\\\
k^2&=&\frac{5{e_1}^2}{1-{e_1}^2}\frac{1-h}{h+4{e_1}^2}\ .\end{aligned}$$

The constant $k^2$ is positive because of relation . The value $k^2=1$ defines a limit between two dynamical regimes. If $k^2<1$, or equivalently if $h>{e_1}^2$, $\dot{\phi}$ never vanishes and the projection of ${\mathbf{w}}$ on the orbital plane of the binary moves along the full ellipse . In this case, ${\mathbf{w}}$ precesses around the angular momentum of the binary, ${\mathbf{k}}$. If $z>0$ the mutual inclination of the two orbits is always less than $90^\circ$ so the orbital motion of the particle is prograde; conversely, if $z<0$ the mutual inclination of the two orbits is always superior to $90^\circ$ so the orbital motion of the particle is retrograde.

If $k^2>1$ (or $h<{e_1}^2$), then $\dot{\phi}$ vanishes for $\phi_0=\pm \arcsin(1/k)$, changes its sign (which is accompanied by a change of sign in the $z$ variable), and the angle $\phi$ librates between $-\phi_0$ and $+\phi_0$. Thus, the projection of ${\mathbf{w}}$ on the orbital plane of the binary is bounded by the unit circle to stay on an arc of ellipse . In this case, ${\mathbf{w}}$ precesses around the direction of perihelion of the binary, so that both the inclination and the node of the particle librate around $\pm 90^\circ$.

In both cases, the period of the motion can be calculated with a simple quadrature using equation :

$$\label{eq:Tpart}
 T=\frac{16}{3n_1}\frac{{\mathcal{M}_{01}}}{{\beta_1}}\left(\frac{a_2}{{a_1}}\right)^{7/2}
 \frac{K(k^2)(1-e_2^2)^2}{\sqrt{(1-{e_1}^2)(h+4{e_1}^2)}}\ ,$$

where $K(k^2)$ is the elliptic integral of the first kind defined by: $$K(k^2) = \left\{\begin{array}{cc}
                  \int_0^{\pi/2}\frac{d\phi}{\sqrt{1-k^2\sin^2\phi}} & \text{if}\ k^2<1\\
                  \int_0^{\phi_0}\frac{d\phi}{\sqrt{1-k^2\sin^2\phi}} & \text{if}\ k^2>1
                 \end{array}\right.\ .$$

The last case where $k^2=1$ (or $h={e_1}^2$) corresponds to the trajectories that separate the previous two types. They link the points $(x,y,z)=(0,\pm 1,0)$, and the associated period is infinite. In the projection on the $(x,y)$ plane, these separatrices form the ellipse which is tangent to the unit circle. Since all trajectories that are inside this ellipse correspond to the precession of ${\mathbf{w}}$ around ${\mathbf{k}}$, the width $\Delta x_{\text{sep}}$ of the separating ellipse in the $(x,y)$ plane gives an indication on the proportion of such trajectories. Using equation and the fact that $h={e_1}^2$ on the separatrix, we get:

$$\Delta x_{\text{sep}} = 2\sqrt{\frac{1-{e_1}^2}{1+4{e_1}^2}}$$

Therefore, when the inner binary is circular, this width is equal to $2$, the full width of the unit circle, and the only possible motion is precession of ${\mathbf{w}}$ around $\pm{\mathbf{k}}$. When the eccentricity of the binary increases, the width of the separatrix decreases to zero, which is a limit case since it can only be reached for a value of the binary’s eccentricity equal to 1. The precession motions of ${\mathbf{w}}$ around $\pm{\mathbf{i}}$ thus become predominant when the eccentricity of the binary grows.

Comparison with numerical studies
---------------------------------

![Energy levels of the Hamiltonian in the $(i_2\cos\Omega_2,i_2\sin\Omega_2)$ plane for values of the eccentricity of the binary ${e_1}=0$ (a), ${e_1}=0.1$ (b), ${e_1}=0.79$ (c), ${e_1}=0.9$ (d).\[fig:Heb\]](FL2009-fig2){width="\columnwidth"}

In [@2009MNRAS.394.1721V], the authors investigate a family of particles at high inclinations around the binary HD 98800 Ba-Bb, which remain stable even under the perturbation of an outer third stellar companion. They isolate a nodal precession imposed by the inner binary as the stabilizing mechanism working against the destabilizing Kozai perturbations of the outer companion. They run simulations of test particles orbiting the binary HD 98800 Ba-Bb using non secular equations. They observe the libration islands around $i_2=\pm90^\circ$ and $\Omega_2=\pm90^\circ$ that we discussed in the previous section. As they show their results in the $(i_2\cos\Omega_2,i_2\sin\Omega_2)$ plane, we plotted the energy levels of the outer restricted Hamiltonian using these same coordinates for an easier comparison. Figure \[fig:Heb\] shows these levels for different values of the eccentricity. The c. panel in particular uses the same value for the eccentricity of the binary (${e_1}=0.79$) as figures 4 and 5 of [@2009MNRAS.394.1721V].

Verrier and Evans notice no apparent structure in the dynamics of the $(e_2,\omega_2)$ couple apart from the circulation of the perihelion. This is in agreement with the fact that the particle’s eccentricity is constant at the quadrupolar approximation.

The authors also suggest that the projection of the angular momentum of test particles along the line of apses of the binary may be an integral of motion. From the results of the previous section, it is straightforward to see that the projection $x$ of the angular momentum of test particles along the line of apses of the binary is not constant. It varies with an amplitude that decreases to $0$ when the inclination of the particle approaches $\pm 90^\circ$, which can be misleading when looking at numerical results for highly inclined particles. However, the norm of the angular momentum of the test particles is an integral of the secular motion.

![\[fig:period\] Dependence of the period with respect to the eccentricity of the binary, in normalized units. The full line corresponds to the calculated period, while the dashed line corresponds to the power-law fit given by @2009MNRAS.394.1721V. We used a least squares method to fit the relative position of the two curves.](FL2009-fig3){width="\columnwidth"}

Finally, the authors give a power-law fit of the period of the libration of the node with respect to three parameters: the eccentricity of the binary, the ratio of the semimajor axes $a_2/a_1$, and the mass ratio of the binary, $\delta$. Their power-law is fitted using particles with fixed inclinations ($85^\circ$). They give in their equation (5): $$T\propto e_1^{-1.1}\delta^{-0.8}\left(\frac{a_2}{a_1}\right)^{3.37}\ .$$

By rewriting the mass dependences of equation , we get the following analytical dependence with respect to the mass ratio and the semi-major axis of the binary: $$T\propto\left(\frac{a_2}{a_1}\right)^{3.5}(\delta(1-\delta))^{-1}\ .$$

These two exponents compare very well with the fitted power law, in spite of the differences between the two models. The dependency with respect to $e_1$ is rather complex in equation , and it is best compared in figure \[fig:period\].

The grid of initial conditions for the particles in [@2009MNRAS.394.1721V] extends however from 3 to 10 AU for a binary separation of 1 AU, so the quadrupolar approximation may not be sufficient to fully describe the motion of the particles with the lowest semi major axes. In particular, Verrier and Evans state that some low inclination particles show large eccentricity variations and even instability. This could be due to a low initial semi-major axis and to resonances that are eliminated in our secular model by the averaging over the mean anomalies.

Problem of three massive bodies {#se:ftbp}
===============================

As we already stated, the quadrupolar secular three-body problem is still integrable when all the bodies have positive masses. As such, it is possible to show how the outer restricted problem we discussed in the previous section relates to the general case, and to the inner restricted case studied by @1962AJ.....67..591K and . We will first express the Hamiltonian of the secular quadrupolar problem using the same vectorial method as in the previous section in order to focus on the relative movements of the orbits. In their studies of the secular quadrupolar problem, @1976CeMec..13..471L and @1994CeMDA..58..245F have shown that this problem depends on two parameters. We will then point out which regions of parameter space are topologically equivalent to the outer restricted case, and which regions correspond to the inner restricted case, in order to show the continuity that exists between both situations.

Hamiltonian
-----------

Let us consider three masses $m_0$, $m_1$ and $m_2$, this time with $m_2\neq0$. We note the barycentric coordinates and impulsions $({\mathbf{u}}_i,{\mathbf{\tilde{u}}}_i)_{i=0,1,2}$. As in the previous section, we suppose that the two bodies of indices 0 and 1 form a binary and that the distance of the third body to this binary is much larger than the separation of the binary. We still note $\delta =
m_0/(m_0+m_1)$. We first perform a canonical change of variables to Jacobi coordinates,

$$\begin{aligned}
 {\mathbf{r}}_0 &={\mathbf{u}}_0 &\tilde{{\mathbf{r}}}_0 &= \tilde{{\mathbf{u}}}_0+\tilde{{\mathbf{u}}}_1+\tilde{{\mathbf{u}}}_2 ={\mathbf{0}} \\
 {\mathbf{r}}_1 &={\mathbf{u}}_1-{\mathbf{u}}_0 &\tilde{{\mathbf{r}}}_1 &= {\mathbf{\tilde{u}}}_1 + (1-\delta)\tilde{{\mathbf{u}}}_2 \\
 {\mathbf{r}}_2 &={\mathbf{u}}_2-\delta{\mathbf{u}}_0-(1-\delta){\mathbf{u}}_1 & \tilde{{\mathbf{r}}}_2 &= {\mathbf{\tilde{u}}}_2\end{aligned}$$

In these coordinates, the Hamiltonian of the three bodies is [@1989LaskarGoutelas]: $$\begin{gathered}
 H=\left(\frac{{\tilde{\mathbf{r}}}_1^2}{2\beta_1}-\frac{\mu_1\beta_1}{r_1}\right) + \frac{{\tilde{\mathbf{r}}}_2^2}{2\beta_2}\\
   -Gm_2\left(\frac{m_0}{|{\mathbf{r}}_2+(1-\delta){\mathbf{r}}_1|}+\frac{m_1}{|{\mathbf{r}}_2-\delta{\mathbf{r}}_1|}\right),\end{gathered}$$

where $\beta_1^{-1}=m_0^{-1}+m_1^{-1}$, $\beta_2^{-1}=(m_0+m_1)^{-1}+m_2^{-1}$, $\mu_1=G(m_0+m_1)$ and $\mu_2=G(m_0+m_1+m_2)$.

Using the fact that $r_1\ll r_2$, we expand the Hamiltonian to order two in $r_1/r_2$ as in the previous section: $$\begin{gathered}
 H=\left(\frac{{\tilde{\mathbf{r}}}_1^2}{2\beta_1}-\frac{\mu_1\beta_1}{r_1}\right)+
 \left(\frac{{\tilde{\mathbf{r}}}_2^2}{2\beta_2}-\frac{\mu_2\beta_2}{r_2}\right)\\
 -G\frac{\beta_1m_2}{2r_2^3}\left(3\frac{({\mathbf{r}}_1.{\mathbf{r}}_2)^2}{r_2^2}-r_1^2\right)\ .\end{gathered}$$

The first two terms are Keplerian energies and are equal respectively to $-\mu_1\beta_1/2a_1$ and $-\mu_2\beta_2/2a_2$, where $a_1$ and $a_2$ are the semi major axes of the inner and the outer body in our system of coordinates.

We now average over the two mean anomalies $M_1$ and $M_2$ in order to get the secular part of the Hamiltonian. We will define 4 unit vectors which are analogous to the 4 vectors we used in the first section: $({\mathbf{i}}_1,{\mathbf{j}}_1,{\mathbf{k}}_1)$ are tied to the orbit of the inner binary; ${\mathbf{i}}_1$ points in the direction of the perihelion of the inner binary, ${\mathbf{k}}_1$ points in the direction of its angular momentum, and ${\mathbf{j}}_1={\mathbf{k}}_1\wedge{\mathbf{i}}_1$. The last vector ${\mathbf{k}}_2$ is colinear to the angular momentum of the outer body. In this section, the vectors tied to the orbit of the inner binary will no longer have fixed directions.

Using the same averaging formulae as in the previous section and using the fact that $({\mathbf{i}}_1.{\mathbf{k}}_2)^2+
({\mathbf{j}}_1.{\mathbf{k}}_2)^2+({\mathbf{k}}_1.{\mathbf{k}}_2)^2={\mathbf{k}}_2^2=1$, we can write:

$$\begin{gathered}
\label{eq:tmpavg2}
  {\left\langle H \right\rangle}_{M_1,M_2}=-\frac{\mu_1\beta_1}{2a_1}-\frac{\mu_2\beta_2}{2a_2}
  -\frac{3}{8}\frac{Gm_2\beta_1}{(1-e_2^2)^{3/2}}\frac{a_1^2}{a_2^3}\\
  \times\left[-\frac{1}{3}+2e_1^2+
  (1-e_1^2){\mathopen ( {\mathbf{k}}_1 . {\mathbf{k}}_2 \mathclose )}^2-5e_1^2{\mathopen ( {\mathbf{i}}_1 . {\mathbf{k}}_2 \mathclose )}^2\right]\ .\end{gathered}$$

After averaging over the two mean anomalies, the semi-major axes are constant. There are thus four degrees of freedom in the system, associated to the two eccentricities and the two inclinations. As we explained in the previous section, the argument of perihelion of the outer body does not appear in the quadrupolar expansion, and thus the norm of the angular momentum of the outer body $G_2=\beta_2\sqrt{\mu_2a_2(1-e_2^2)}$, is constant. This implies that its eccentricity $e_2$ is constant. Using the reduction of the nodes would leave only one degree of freedom in the reduced Hamiltonian, associated to the couple $(e_1,\omega_1)$. The full reduction of the Hamiltonian and its expression in elliptical variables is the approach that has been used widely, since it yields a very similar Hamiltonian function as in the inner restricted problem [@1969CeMec...1..200H; @1976CeMec..13..471L; @1994CeMDA..58..245F].

We want however to look at the motion of the nodes, or equivalently the motion of the vector ${\mathbf{k}}_2$ in the moving frame $({\mathbf{i}}_1,{\mathbf{j}}_1,{\mathbf{k}}_1)$ of the orbit of the second body.

In order to easily compute the equations of motion, we introduce two vectors associated to the orbit of the binary that are colinear to ${\mathbf{i}}_1$ and ${\mathbf{k}}_1$, and include in their norm the eccentricity of the binary, as in [@2009AJ....137.3706T]: $$\begin{aligned}
 {\mathcal{K}}_1 &= \sqrt{1-e_1^2}{\mathbf{k}}_1\ , &{\mathcal{I}}_1 &=e_1 {\mathbf{i}}_1\ .\end{aligned}$$

If we drop all the constant terms in equation \[eq:tmpavg2\] and use the above vectors, we get:

$${\left\langle H \right\rangle}_{M_1,M_2}=-\frac{\alpha'G_2}{2}\left[2
 {\mathcal{I}}_1^2+({\mathcal{K}}_1.{\mathbf{k}}_2)^2
                  -5({\mathcal{I}}_1.{\mathbf{k}}_2)^2\right]\ ,$$

where $$\alpha' = \frac{3}{4}n_1\left(\frac{a_1}{a_2}\right)^{7/2}\frac{\beta_1}{{\mathcal{M}_{01}}}\frac{1}{(1-e_2^2)^2}
        \sqrt{1+\frac{m_2}{{\mathcal{M}_{01}}}}\ ,$$ and ${\mathcal{M}_{01}},n_1$ are defined as in section \[se:rtbp\].

Equations of motion {#equations-of-motion}
-------------------

The components of ${\mathcal{K}}_1,{\mathcal{I}}_1$ and ${\mathbf{k}}_2$ have the following Poisson brackets[^4] [@2005BorisovMamaev; @2006Icar..185..312B]:

$$\begin{aligned}
 \{{\mathcal{K}}_{1i} ,{\mathcal{K}}_{1j} \} &= -\frac{\epsilon_{ijk}}{\Lambda_1} {\mathcal{K}}_{1k}\,, & 
           \{k_{2i} ,k_{2j} \} &= -\frac{\epsilon_{ijk}}{G_2} k_{2k}\,, \\
 \{{\mathcal{I}}_{1i} ,{\mathcal{I}}_{1j} \} &= -\frac{\epsilon_{ijk}}{\Lambda_1} {\mathcal{K}}_{1k}\,, & 
           \{{\mathcal{K}}_{1i} ,{\mathcal{I}}_{1j} \} &= -\frac{\epsilon_{ijk}}{\Lambda_1} {\mathcal{I}}_{1k}\,,\end{aligned}$$

where $\Lambda_1=\beta_1\sqrt{\mu_1a_1}$ and $\epsilon_{ijk}$ is the Levi-Civita symbol[^5].

The equations of motion for the three vectors are thus: $$\begin{aligned}
 \dot{{\mathcal{K}}}_1 & = & - \frac{1}{\Lambda_1} \left( {\mathcal{K}}_1 \wedge \nabla_{{\mathcal{K}}_1} H 
                      + {\mathcal{I}}_1 \wedge \nabla_{{\mathcal{I}}_1} H\right)\ , \\
 \dot{{\mathcal{I}}}_1 & = & - \frac{1}{\Lambda_1} \left( {\mathcal{I}}_1 \wedge \nabla_{{\mathcal{K}}_1} H 
                      + {\mathcal{K}}_1 \wedge \nabla_{{\mathcal{I}}_1} H\right)\ , \\
 \dot{{\mathbf{k}}}_2 & = & - \frac{1}{G_2} {\mathbf{k}}_2 \wedge \nabla_{{\mathbf{k}}_2} H\ .\end{aligned}$$

In order to look at the motion of the vector ${\mathbf{k}}_2$ in the moving frame $({\mathbf{i}}_1,{\mathbf{j}}_1,{\mathbf{k}}_1)$ of the orbit of the second body, we use as @2006Icar..185..312B the above system to derive equations for $x=({\mathbf{k}}_2.{\mathbf{i}}_1)$, $y=({\mathbf{k}}_2.{\mathbf{j}}_1)$, $z=({\mathbf{k}}_2.{\mathbf{k}}_1)$ and $e_1$. Indeed, $x=({\mathbf{k}}_2.{\mathcal{I}}_1)/|{\mathcal{I}}_1|$, $z=({\mathbf{k}}_2.{\mathcal{K}}_1)/|{\mathcal{K}}_1|$, $e_1=|{\mathcal{I}}_1|$, and $y$ is obtained using the identity $x^2+y^2+z^2=1$:

$$\begin{aligned}
 \dot{x} & =  \alpha'(1-e_1^2)yz + \alpha'\frac{G_2}{\Lambda_1}\sqrt{1-e_1^2}y(2-5x^2) 
 \label{eq:fullmotion1}\\
 \dot{y} & =  -\alpha'(1+4e_1^2)xz  \nonumber \\
 &\phantom{=\dots}-\alpha'\frac{G_2}{\Lambda_1}\frac{x}{\sqrt{1-e_1^2}}[(1-e_1^2)(2-5x^2)+5e_1^2z^2] \\
 \dot{z} & =  5\alpha'e_1^2xy +\alpha'\frac{G_2}{\Lambda_1}\frac{5e_1^2}{\sqrt{1-e_1^2}}xyz \\
 \dot{e}_1 &=  \alpha'\frac{G_2}{\Lambda_1}5e_1\sqrt{1-e_1^2}xy \label{eq:fullmotion4}\end{aligned}$$

The equations for $x,y$ and $z$ contain two terms: the first one is identical to the outer restricted system, and the second one includes the motion of the reference frame $({\mathbf{i}}_1,{\mathbf{j}}_1,{\mathbf{k}}_1)$ induced by the interaction with the third body. Note that when $\Lambda_1$ is very large compared to $G_2$ so that we can assume that $G_2/\Lambda_1$ is equal to zero, which corresponds to the case where $m_2\ll m_0\text{ and }m_1$, the above system is identical to the outer restricted system –.

The conservation of the total angular momentum ${\mathbf{C}}={\mathbf{G}}_1+{\mathbf{G}}_2$ introduces the two main parameters of the problem. Indeed, $$\Lambda_1^2(1-e_1^2) + G_2^2 + 2 \Lambda_1\sqrt{1-e_1^2} G_2 z = C^2\ .$$

We note $\gamma=C/\Lambda_1$, $\gamma_2=G_2/\Lambda_1$. The above expression of the norm of the total angular momentum can be rewritten as a second degree equation giving $\sqrt{1-e_1^2}$ as a function of $z$ using the two parameters $\gamma$ and $\gamma_2$: $$\label{eq:redangmom}
 (1-e_1)^2 + 2\gamma_2z\sqrt{1-e_1^2} + \gamma_2^2-\gamma^2 = 0\ .$$

The Hamiltonian can then be rewritten as: $${\left\langle H \right\rangle}=-\frac{1}{2}\alpha'\Lambda_1\gamma_2[z^2+e_1^2(2-z^2-5x^2)]\ .$$

The inequalities $-1\leqslant z\leqslant 1$ and $0\leqslant e_1 < 1$ give the boundaries of the parameter space and the range of possible values of $e_1$ for any given couple of parameters[^6] $(\gamma,\gamma_2)$: $$\begin{aligned}
    |\gamma-\gamma_2|&\leqslant 1\ , \label{eq:paramspace}\\
    |\gamma-\gamma_2|&\leqslant \sqrt{1-e_1^2} \leqslant \min[\gamma+\gamma_2,1] \label{eq:allowedeta}\ . \end{aligned}$$

With these notations, the outer restricted problem of section \[se:rtbp\] corresponds to the limit where $\gamma_2=0$, and in this case $e_1=\sqrt{1-\gamma^2}$ is constant as we saw. Note that when $\gamma_2>\gamma$, we have $G_2>C$, so this part of the parameter space contains only retrograde motions. Our aim in this paper is to show the continuity between the outer restricted case we studied in section \[se:rtbp\], and the inner restricted case that was investigated by @1962AJ.....67..591K and . Both these problems lie in the region of parameter space where $\gamma>\gamma_2$ so we will restrict our study to this case[^7].

In our case where $\gamma>\gamma_2$, there is only one acceptable root to [equation ]{}, which is: $$\sqrt{1-e^2_1}=-\gamma_2z+\sqrt{(\gamma_2z)^2+\gamma^2-\gamma_2^2}\ .$$

This relation implies that $e_1$ is a growing function of $z$. Note that $z=\cos i_2$, where $i_2$ is the inclination of the outer body in the reference frame of the inner binary. As such, coplanar prograde motions ($z=1$) will always occur for the maximal value of the eccentricity of the inner binary: $$e_{1,\max}=\sqrt{1-(\gamma-\gamma_2)^2}\ .$$

Conversely, low eccentricities for the binary will be associated to lower values of $z$, and thus higher inclinations. Relation implies that the inner binary can only have a circular motion if $\gamma+\gamma_2\geqslant 1$. In this case, coplanar retrograde motion ($z=-1$) is not allowed, and the lowest value of $z$ is: $$\label{eq:z0}
   z_0 = \cos i_{2,\max} = \frac{\gamma^2-\gamma_2^2-1}{2\gamma_2}\ .$$

When $\gamma+\gamma_2<1$ however, coplanar retrograde motion ($z=-1$) is possible and the associated value of the eccentricity of the binary is: $$e_{1,\min}=\sqrt{1-(\gamma+\gamma_2)^2}\ .$$

![\[fig:paramspace\] Parameter Space. The dark gray areas are excluded by equation , the light gray area correspond to the part of parameter space corresponding to $\gamma_2\geqslant\gamma$ which we do not study. The dotted line $\gamma+\gamma_2=1$ separates the zone where there can be coplanar retrograde motion associated to a minimum eccentricity for the inner binary that is strictly higher than $0$ (below the dotted line) and the zone where the inner binary can be circular but the inclination is bounded (see section \[sse:eqcirc\]). In zone O, the problem is topologically equivalent to the outer restricted problem. In zones I and I’ it is topologically equivalent to the inner restricted problem, with zone I being equivalent to situations above the critical inclination and zone I’ being equivalent to situations under the critical inclination. The letters a–f correspond to the values of the parameters used to plot the corresponding panels in figures \[fig:IOPlane\] and \[fig:Spheres\].](FL2009-fig4){width="\columnwidth"}

Dynamical regimes {#se:dynreg}
-----------------

![\[fig:IOPlane\] Trajectories in the $(i_2\cos\Omega_2,i_2\sin\Omega_2)$ plane for different values of the parameters. [See section \[se:dynreg\] for a detailed discussion and appendix \[se:fixp\] for calculations.]{} a: $(\gamma,\gamma_2)=(0.8,0)$ outer restricted case with $e_1=0.6$; b: $(\gamma,\gamma_2)=(0.8,0.25)$; c: $\gamma=0.8,\,
\gamma_2^2=(1/3)(1-\gamma^2)$; d: $(\gamma,\gamma_2)=(0.8,0.4)$; e: $(\gamma,\gamma_2)=(1.08,0.4)$; f: $(\gamma,\gamma_2)=(1.28,0.4)$.](FL2009-fig5){width="\columnwidth"}

![\[fig:Spheres\] Trajectories on the unit angular momentum sphere for different values of the parameters. [See section \[se:dynreg\] for a detailed discussion and appendix \[se:fixp\] for calculations.]{} a: $(\gamma,\gamma_2)=(0.8,0)$ outer restricted case with $e_1=0.6$; b: $(\gamma,\gamma_2)=(0.8,0.25)$; c: $\gamma=0.8,\,
\gamma_2^2=(1/3)(1-\gamma^2)$; d: $(\gamma,\gamma_2)=(0.8,0.4)$; e: $(\gamma,\gamma_2)=(1.08,0.4)$; f: $(\gamma,\gamma_2)=(1.28,0.4)$.](FL2009-fig6){width="\columnwidth"}

In appendix \[se:fixp\], we briefly derive in the framework of the present study the fixed points of the system and the boundaries of the dynamical regimes in parameter space that are given in [@1994CeMDA..58..245F]. [The fixed points are named as follow: the north pole is called $N$, and the south pole $S$; linearly stable fixed points are named $E$, as elliptic, and linearly unstable points are named $H$, as hyperbolic; finally, signs are placed as indices to refer to the symmetry of the problem with respect to the two planes $x=0$ and $y=0$.]{} There are three dynamical regimes in the region of parameter space we study.

In region O of figure \[fig:paramspace\], the parameter $\gamma_2=G_2/\Lambda_1$ is small (less than 1/2). This is the case in particular when the mass ratio $m_2/m_1$ is small. [Moreover, $\gamma^2+3\gamma_2^2<1$.]{} The phase space is topologically equivalent to the outer restricted problem of section \[se:rtbp\]. The north pole, which corresponds to coplanar prograde motion with maximal eccentricity for the binary is linearly stable. There are two additional stable fixed points $\text{E}_{\pm}$ in the plane $y=0$ [(see section \[sse:eqy0\])]{}. They belong to the same family as the fixed points $y=z=0$, $x=\pm 1$ of the outer restricted problem that are responsible for the stable high [inclination]{} orbits observed by @2009MNRAS.394.1721V. When $\gamma+\gamma_2\leqslant1$, the south pole which corresponds to coplanar retrograde motion with minimal eccentricity for the binary, is also [linearly]{} stable. [The a panels of figures \[fig:IOPlane\] and \[fig:Spheres\] provide a visualisation of the topology of this case.[^8]]{} When $\gamma+\gamma_2>1$, the south pole is no longer accessible as stated in the previous section. It [is however replaced by a stable trajectory]{} at a [maximal]{} inclination given by equation , as can be seen on panel b [of figures \[fig:IOPlane\] and \[fig:Spheres\]]{}. This trajectory corresponds to a circular inner binary [(see section \[sse:eqcirc\])]{}. Finally, there are two unstable points $\text{H}_{\pm}$ in the $x=0$ plane that belong to the same family as the unstable points of the outer restricted problem $x=z=0$, $y=\pm 1$ [(see section \[sse:eqx0\])]{}.

[Panels c of figures \[fig:IOPlane\] and \[fig:Spheres\] show the limiting case between regions O and I. On this boundary, $\gamma^2+3\gamma_2^2=1$. The two unstable points $H_{\pm}$ are now located on the boundary of the accessible part of the sphere.]{}

[Regions I and I’ of figure \[fig:paramspace\] are both in the part of the parameter space defined by $\gamma^2+3\gamma_2^2>1$. In this zone, the problem becomes topologically equivalent to the inner restricted problem studied by and @1962AJ.....67..591K. In the inner restricted case, there is a critical value of the inclination ($\cos i_2=\sqrt{3/5}$) under which a circular inner binary is always linearly stable, and above which a circular inner binary is always linearly unstable, giving rise to Kozai cycles.]{}

In region I of figure \[fig:paramspace\], the dynamical regime is topologically equivalent to the inner restricted problem in the case where the inclination is superior to the critical value. The limit trajectory $z=z_0$ which corresponds to a circular inner binary becomes [linearly]{} unstable. However, the north pole and the two fixed points $\text{E}_{\pm}$ are still stable. In the inner restricted phase space, when the inclination is superior to the critical value, there are two possible behaviours for the periastron of the inner particle: it can either circulate, or librate around $\pm 90^{\circ}$. In our representation, the circulation case corresponds to trajectories around the north pole, and the libration islands correspond to the two fixed points $\text{E}_{\pm}$. This is shown in panels d and e in figures \[fig:IOPlane\] and \[fig:Spheres\].

In region I’ of figure \[fig:paramspace\], the dynamical regime is topologically equivalent to the inner restricted problem in the case where the inclination is inferior to the critical value. Only one stable fixed point remains, on the north pole of the sphere, associated to prograde coplanar motion. This is shown in panel f in figures \[fig:IOPlane\] and \[fig:Spheres\].

In both regions I and I’, the parameter $\gamma_2=G_2/\Lambda_1$ can take higher values. This is in particular true when the mass ratio $m_2/m_1$ increases.

The curve between regions I and I’ is linked to the critical inclination that is defined in the inner restricted case. Indeed, along that curve, given by equation \[eq:critlim\], we have the following limits when $\gamma \rightarrow \infty$: $$\begin{aligned}
    \frac{\gamma_2}{\gamma} &\rightarrow 1\, , & \gamma_2-\gamma \rightarrow -\sqrt{\frac{3}{5}}\ . \end{aligned}$$

When $G_2$ is very large compared to $G_1$, we can make the following first order expansion:

$$\begin{aligned}
 \gamma_2-\gamma &=\frac{G_2-C}{\Lambda_1} \\
  &=\frac{G_2-\sqrt{G_2^2+G_1^2+2({\mathbf{G}}_2.{\mathbf{G}}_1)}}{\Lambda_1} \\
  &\approx \frac{G_2 - G_2(1+ ({\mathbf{G}}_2.{\mathbf{G}}_1)/G_2^2)}{\Lambda_1} \\
  &\approx - \frac{({\mathbf{G}}_2.{\mathbf{G}}_1)}{\Lambda_1G_2} \\
  &\approx -z \sqrt{1-e_1^2} \ .\end{aligned}$$

As such, we see that along the border between regions I and I’, when $\gamma$ and $\gamma_2$ both tend to infinity, we have the relation

$$z \sqrt{1-e_1^2} \approx \sqrt{\frac{3}{5}}\ .$$

Recall that $z=\cos i_2$, where $i_2$ is the inclination of the outer orbit in the reference frame of the inner orbit. Thus, the inclination of the inner orbit relatively to the outer orbit is $i_1=-i_2$, and the above equation becomes:

$$\cos i_1 \sqrt{1-e_1^2} \approx \sqrt{\frac{3}{5}}\ .$$

This relation is precisely the one giving the critical value of the normal component of the angular momentum of the inner body in the inner restricted problem.

Conclusion {#conclusion .unnumbered}
==========

We first studied the case of a massless particle orbiting a binary at a long distance, and, in the secular and quadrupolar approximations, gave a full analytical description of the motion along with the expression of the period of the secular motion. When the inner binary is circular, only nodal precession takes place. However, when the binary is elliptic, libration islands appear at high inclinations, and these islands grow bigger when the eccentricity of the binary rises. @2008MNRAS.390.1377V [@2009MNRAS.394.1721V] observe a similar nodal libration in their study of the stability of particle populations in the quadruple stellar system HD 98800, and we showed that the analytical framework that we derived for the outer restricted problem is well suited to explain the results of Verrier and Evans.

The quadrupolar secular three-body problem is still integrable when all the bodies have positive masses [@1969CeMec...1..200H; @1976CeMec..13..471L; @1994CeMDA..58..245F]. Using a vectorial formalism as [@2006Icar..185..312B; @2009Icar..201..750B; @2009AJ....137.3706T], we looked at this problem from the point of view of the outer restricted case. We showed how the outer restricted problem relates to the general case, and to the inner restricted case studied by @1962AJ.....67..591K and : when the mass of the outer body is small enough compared to the mass of the inner body, the general case behaves similarly to the outer restricted problem. When the mass of the outer body increases enough, the general case behaves like the inner restricted problem. We gave an expression of the boundary between these two regimes.

The outer restricted problem and its generalization to the non restricted case provide an interesting starting point in the study of circumbinary planetary systems, such as the one discovered recently around the eclipsing sdB+M system HW Virginis [@2009AJ....137.3181L]. In this system, the inner binary is very tight with a period of 2.8 hr, while the proposed planetary companions have periods of 9.1 yr and 15.8 yr, so the quadrupolar expansion is fully justified. Another field of application of the outer restricted problem is the study of the motion of stars orbiting around binary black holes [@2008AJ....135.2398M; @2009ApJ...692.1075G; @2009ApJ...693L..35M].

Fixed points and bifurcations {#se:fixp}
=============================

The fixed points and boundaries presented in section \[se:dynreg\] have already been studied by @1976CeMec..13..471L and @1994CeMDA..58..245F. We briefly present here their derivation in the framework of the present formalism. With the notations of section \[se:ftbp\], we will limit ourselves to $\gamma>\gamma_2$.

Poles of the sphere, $x=y=0$ {#sse:eqpol}
----------------------------

This case corresponds to case 1 in section 5 of [@1994CeMDA..58..245F]. Note that their sphere is constructed using the eccentricity and perihelion of the [inner binary]{}, and is thus different from our angular momentum sphere.

For all values of the parameters in the domain we study, the north pole $z=1$, which corresponds to coplanar prograde motion, is a fixed point of the system. The associated eccentricity of the inner binary is: $$e_{1,\max}=\sqrt{1-(\gamma-\gamma_2)^2}\ .$$ It is the maximal value of the eccentricity. This fixed point is always linearly stable. It is noted N in figures \[fig:IOPlane\] and \[fig:Spheres\]. Figure \[fig:IOPlane\] shows the lines of equal energy in the $(i_2\cos\Omega_2,i_2\sin\Omega_2)$ plane, and figure \[fig:Spheres\] shows these lines plotted on the sphere of unit angular momentum of the outer body ${\mathbf{k}}_2^2=1$.

When $\gamma+\gamma_2<1$ (under the dotted line in figure \[fig:paramspace\]), the south pole $z=-1$ (noted S in the following figures), which corresponds to coplanar retrograde motion, is also a linearly stable fixed point of the system. The eccentricity of the inner binary is minimal and equal to: $$e_{1,\min}=\sqrt{1-(\gamma+\gamma_2)^2}\ .$$

Note that in this region of parameter space, the inner binary cannot be circular.

When $\gamma+\gamma_2\geqslant1$ (above the dotted line in figure \[fig:paramspace\]), the minimal eccentricity [of]{} the binary is $0$ as deduced from . The south pole $z=-1$ does not correspond to a real value of the eccentricity in this case. This limit however is not a bifurcation strictly speaking. When crossing it, the stable south pole of the sphere [is replaced by a stable trajectory at maximal inclination]{}.

Circular Trajectories for the inner binary {#sse:eqcirc}
------------------------------------------

In the region of parameter space where circular trajectories exist for the binary (above the dotted line in figure \[fig:paramspace\]), the value of $z$ which corresponds to such trajectories is minimal and equal to: $$z_0= \frac{\gamma^2-\gamma_2^2-1}{2\gamma_2}\ .$$

The equations of motion on the small circle of the sphere $z=z_0$ are: $$\begin{aligned}
 \dot{x}_{\phantom{1}} & = & \phantom{-} \alpha' y(z_0 + \gamma_2 (2-5x^2))\ , \label{eq:circmotionx}\\
 \dot{y}_{\phantom{1}} & = & - \alpha'x(z_0 +\gamma_2 (2-5x^2))\ , \label{eq:circmotiony}\\
 \dot{z}_{\phantom{1}} & = & \phantom{-} 0\ , \\
 \dot{e}_1 &= &\phantom{-}  0\ .\end{aligned}$$

The right hand sides of equations and vanish for a certain value of $x$ equal to: $$x_0^2 = \frac{\gamma^2+3\gamma_2^2-1}{10\gamma_2^2}\ .$$

The curve $\gamma^2+3\gamma_2^2=1$ separates in figure \[fig:paramspace\] the regions noted O and I. We can distinguish three cases:

1.  $\gamma^2+3\gamma_2^2<1$. In region O of figure \[fig:paramspace\], $x_0^2<0$ so there is no fixed point on the circle $z=z_0$. As such, this circle is a trajectory for which the inner binary is circular and the outer orbit precesses at a fixed inclination given by $i_{2,\max}=\text{ArcCos}\,z_0$. Moreover, this trajectory is linearly stable.

2.  $\gamma^2+3\gamma_2^2=1$. There are two fixed points on the circle $z=z_0$ at the coordinates $(x=0,y=\pm\sqrt{1-z_0^2})$.

3.  $\gamma^2+3\gamma_2^2> 1$. In this case, we must also check that $y_0^2=1-x_0^2-z_0^2\geqslant 0$. The limit case where there is equality yields: $$5 \gamma_2^4 - (4 + 10 \gamma^2)\gamma_2^2 +  (5 \gamma^4  - 8 \gamma^2 + 3) = 0 \ .$$

This boundary limits the regions I and I’ in figure \[fig:paramspace\]. When solving the above equation for $\gamma_2^2$ and selecting only the relevant solution satisfying $\gamma>\gamma_2\, , \gamma+\gamma_2\geqslant
1$, we obtain a solution that corresponds to equation 44 in section 5.2 of [@1994CeMDA..58..245F] and that can be written using our notations as: $$\label{eq:critlim}
 \gamma_2^2 = \frac{2+5\gamma^2-\sqrt{60\gamma^2-11}}{5}\, , \quad \gamma+\gamma_2\geqslant 1\ .$$

In region I, $y_0^2>0$ so there are four fixed points on the circle $z=z_0$, at the coordinates $(\pm x_0, \pm
y_0)$. They are noted $\text{H}_{\pm\pm}$, in panels d and e of figures \[fig:IOPlane\] and \[fig:Spheres\]. Moreover, the trajectories that correspond to circular binaries are unstable in this zone. In region I’ however, $y_0^2<0$ so we are again in a region of parameter space where there are no fixed points on the circle $z=z_0$, and the trajectories associated to circular binaries are again stable.

The $x=0$ plane {#sse:eqx0}
---------------

When $x=0$, the only non trivial equation remaining in system – is $\dot{x}=0$. Looking for a fixed point different from $x=y=0$, we have to solve $$z\sqrt{1-e_1^2}+2\gamma_2=0\ ,$$

which after using relation yields:

$$\begin{aligned}
 e_1 & = &\sqrt{1-\gamma^2-3\gamma_2^2}\ ,\\
x & = & 0 \ , \\
y & = & \pm\frac{\sqrt{\gamma^2-\gamma_2^2}}{\sqrt{\gamma^2+3\gamma_2^2}}\ ,\\
z & = & -\frac{2\gamma_2}{\sqrt{\gamma^2+3\gamma_2^2}}\ .\end{aligned}$$

We thus have two symmetric fixed points in the plane $x=0$. They are noted $\text{H}_{\pm}$ in figures \[fig:IOPlane\] and \[fig:Spheres\]. For these fixed points to exist, the associated eccentricity must be a real number. As such, their domain of existence is the region noted O in figure \[fig:paramspace\]. This is case 2.1 in section 5 of [@1994CeMDA..58..245F].

These two fixed points are linearly unstable in their domain of existence. Note that in the outer restricted problem ($\gamma_2=0$), these fixed points are simply $x=z=0, y=\pm1$.

The $y=0$ plane {#sse:eqy0}
---------------

When $y=0$, the only non trivial equation we must solve is $\dot{y}=0$. Here again, we look for another fixed point than $x=y=0$, thus we have to solve:

$$(1+4e_1^2)z +\frac{\gamma_2}{\sqrt{1-e_1^2}}[(1-e_1^2)(2-5x^2)+5e_1^2z^2]=0\ .$$

Substituting $1-z^2$ in place of $x^2$ and then $\sqrt{1-e_1^2}$ in place of $z$ using , we get: $$(1-e_1^2)^3 - \left(\gamma^2+\frac{1}{2}\gamma_2^2+\frac{5}{8}\right)(1-e_1^2)^2+\frac{5}{8}(\gamma^2-\gamma_2^2)^2=0$$

This equation is the same as equation number 40 in [@1994CeMDA..58..245F]. In our region of parameter space, there is at most one root which satisfies to the constraint . The curve separating the zone where there is one solution and the zone where there is no solution corresponds to the [case where the limit value $e_1=0$ is a solution]{}, and coincides with the boundary between regions I and I’ in figure \[fig:paramspace\] which is given by equation .

When there is a solution, the value of $e_1$ can be translated into a value of $z$ using . Since $y=0$, we get two values of $x=\pm \sqrt{1-z^2}$, and there are thus two symmetric fixed points on the sphere, which are both linearly stable. They are noted $\text{E}_{\pm}$ in figures \[fig:IOPlane\] and \[fig:Spheres\]. When $\gamma_2=0$, these fixed points become simply $y=z=0$, $x=\pm 1$, which are responsible of the stable orbits at high inclination as discussed in the previous sections.

[^1]: E-mail: farago@imcce.fr; laskar@imcce.fr

[^2]: We will from now write ${\left\langle H \right\rangle}$ for the averaged Hamiltonian, omitting the subscripts $M_1,M_2$.

[^3]: [The next non-zero term of the Hamiltonian which is the fourth order in $(a_1/a_2)$ plays an important part in the circular case as has been discussed in detail by @2006CeMDA..95...81P.]{}

[^4]: We use the following convention, where $p_i$ are momenta and $q_i$ positions: $\{f,g\}=\sum_i \left(\frac{\partial
    f}{\partial p_i}\frac{\partial g}{\partial q_i} -\frac{\partial f}{\partial q_i}\frac{\partial g}{\partial
    p_i}\right)$.

[^5]: $\epsilon_{ijk}=+1$ if $(i,j,k)$ is an even permutation of $(1,2,3)$, $\epsilon_{ijk}=-1$ is the permutation is odd, and $\epsilon_{ijk}=0$ in all other cases.

[^6]: The left part of the second inequality is strict if $\gamma=\gamma_2$.

[^7]: The other half of the parameter space ($\gamma\leqslant\gamma_2$) corresponds to retrograde motions which are of less physical interest and much more technical to study using our approach, in particular because equation does not have a unique solution in this case. The interested reader will find a complete discussion of this case in [@1976CeMec..13..471L; @1994CeMDA..58..245F].

[^8]: [Note that the south pole in figure \[fig:IOPlane\] a corresponds to the out-most trajectory; this is an artifact of the coordinate map $(i_2\cos\Omega_2,i_2\sin\Omega_2)$ which sends the south pole of the sphere onto the circle $i_2=180^{\circ}$.]{}
