2/1 resonant periodic orbits in three dimensional planetary systems

Assume a planetary system consisting of a star $S$ and two planets $P_{1}$ and $P_{2}$ with masses $m_{0}$, $m_{1}$ and $m_{2}$, respectively, moving in the space under their mutual gravitational attraction. We consider these bodies as point masses and let $O$ be their center of mass. The three bodies are isolated, thus their center of mass is fixed with respect to any inertial system of reference and its angular momentum vector, ${\bf L}$ , is constant. We now introduce a particular inertial coordinate system, $OXYZ$, such that its origin coincides with their center of mass $O$ and its $Z$-axis is parallel to ${\bf L}$. The system is described by six degrees of freedom defined for instance, by the position vectors of the two planets.

We can reduce the number of degrees of freedom by introducing a suitable rotating frame of reference. Following Michalodimitrakis (1979a), we introduce the rotating frame $Gxyz$, such that:

1. 1.

   Its origin coincides with the center of mass $G$ of the bodies $S$ and $P_{1}$.

2. 2.

   Its $z$-axis is always parallel to the $Z$-axis.

3. 3.

   $S$ and $P_{1}$ move always on $xz$-plane.

We define $x_{i}$, $y_{i}$, $z_{i}$ $(i=1,2,3)$ as the Cartesian coordinates of the bodies with respect to $Gxyz$ and an angle $\theta$ as the angle between $X$ and $x$ axes assuming always that $\theta(0)=0$. According to the definition of $Gxyz$, it always holds that $y_{1}=0$. Considering the normalized gravitational constant $G=1$ and setting $m=m_{0}+m_{1}+m_{2}$ (thus $2\pi$ time units (t.u.) correspond to one year), the Lagrangian of the system in the rotating frame of reference is given by Michalodimitrakis (1979a):

$$\begin{array}[]{l}\mathfrak{L}=\frac{\displaystyle 1}{\displaystyle 2}(m_{0}+m_{1})[a(\dot{x}_{1}^{2}+\dot{z}_{1}^{2}+x_{1}^{2}\dot{\theta}^{2})+\displaystyle b[(\dot{x}_{2}^{2}+\dot{y}_{2}^{2}+\dot{z}_{2}^{2})\\[5.69046pt] \quad+\dot{\theta}^{2}(x_{2}^{2}+y_{2}^{2})+2\dot{\theta}(x_{2}\dot{y}_{2}-\dot{x}_{2}y_{2})]]-V,\end{array}$$

(1)

where $V=-\frac{\displaystyle m_{0}m_{1}}{\displaystyle r_{01}}-\frac{\displaystyle m_{0}m_{2}}{\displaystyle r_{02}}-\frac{\displaystyle m_{1}m_{2}}{\displaystyle r_{12}}$, $a=m_{1}/m_{0}$, $b=m_{2}/m$,

$$\begin{array}[]{l}\displaystyle r_{01}^{2}=(\displaystyle 1+\displaystyle a)^{2}(\displaystyle x_{1}^{2}+z_{1}^{2}),\\[8.5359pt] \displaystyle r_{02}^{2}=(\displaystyle ax_{1}+\displaystyle x_{2})^{2}+y_{2}^{2}+(\displaystyle az_{1}+z_{2})^{2},\\[8.5359pt] \displaystyle r_{12}^{2}=(\displaystyle x_{1}-\displaystyle x_{2})^{2}+y_{2}^{2}+(z_{1}-z_{2})^{2}.\end{array}$$

The angle $\theta$ is ignorable and, consequently, the angular momentum $p_{\theta}=\partial{\mathfrak{L}}/\partial{\dot{\theta}}$ is constant and given by

$$p_{\theta}=(m_{0}+m_{1})[ax_{1}^{2}\dot{\theta}+\displaystyle b[\dot{\theta}(x_{2}^{2}+y_{2}^{2})+(x_{2}\dot{y}_{2}-\dot{x}_{2}y_{2})]]=\textnormal{const.}$$

(2)

Solving the Eq. (2) with respect to $\dot{\theta}\$ we find

$$\dot{\theta}=\frac{\frac{p_{\theta}}{m_{0}+m_{1}}-\displaystyle b(x_{2}\dot{y}_{2}-\dot{x}_{2}y_{2})}{\displaystyle ax_{1}^{2}+\displaystyle b(x_{2}^{2}+y_{2}^{2})}.$$

(3)

As far as the components of the angular momentum vector are concerned, due to the choice of the inertial system it always holds:

$$\begin{array}[]{l}L_{X}=(m_{0}+m_{1})[b(y_{2}\dot{z}_{2}-\dot{y}_{2}z_{2})-\dot{\theta}(ax_{1}z_{1}+bx_{2}z_{2})]=0,\\[8.5359pt] L_{Y}=(m_{0}+m_{1})[a(\dot{x}_{1}z_{1}-x_{1}\dot{z}_{1})+b(\dot{x}_{2}z_{2}-x_{2}\dot{z}_{2})-b\dot{\theta}y_{2}z_{2}]=0\end{array}$$

(4)

and

$$L_{Z}=(m_{0}+m_{1})[b(x_{2}\dot{y}_{2}-\dot{x}_{2}y_{2})+\dot{\theta}[ax_{1}^{2}+b(x_{2}^{2}+y_{2}^{2})]],$$

which coincides with integral (2). The constrains defined by Eqs. (4) give:

$$\begin{array}[]{l}z_{1}=\frac{m_{0}m_{2}}{m_{1}m}\left(\frac{y_{2}\dot{z}_{2}-\dot{y}_{2}z_{2}}{x_{1}\dot{\theta}}-\frac{x_{2}z_{2}}{x_{1}}\right),\\[10.0pt] \dot{z}_{1}=\frac{m_{0}m_{2}}{m_{1}m}\frac{1}{x_{1}}(\dot{x}_{2}z_{2}-x_{2}\dot{z}_{2}-\dot{\theta}y_{2}z_{2})+\frac{\dot{x}_{1}}{x_{1}}z_{1}.\end{array}$$

(5)

Thus, for a particular value of the angular momentum, $p_{\theta}$, the position of the system is defined by the coordinate $x_{1}$ of the planet $P_{1}$ and the position $(x_{2},y_{2},z_{2})$ of the planet $P_{2}$. Namely, the system has been reduced to four degrees of freedom. We can derive the equations of motion from the Lagrangian (1):

$$\begin{array}[]{l}\ddot{x}_{1}=-\frac{m_{0}m_{2}(x_{1}-x_{2})}{(m_{0}+m_{1})[(x_{1}-x_{2})^{2}+y_{2}^{2}+(z_{1}-z_{2})^{2}]^{3/2}}-\frac{m_{0}m_{2}(ax_{1}+x_{2})}{(m_{0}+m_{1})[(ax_{1}+x_{2})^{2}+y_{2}^{2}+(az_{1}+z_{2})^{2}]^{3/2}}\\[8.5359pt] \quad-\frac{(m_{0}+m_{1})x_{1}}{[(1+a)^{2}(x_{1}^{2}+z_{1}^{2})]^{3/2}}+x_{1}\dot{\theta}^{2}\end{array}$$

(6)

$$\begin{array}[]{l}\ddot{x}_{2}=\frac{mm_{1}(x_{1}-x_{2})}{(m_{0}+m_{1})[(x_{1}-x_{2})^{2}+y_{2}^{2}+(z_{1}-z_{2})^{2}]^{3/2}}-\frac{mm_{0}(ax_{1}+x_{2})}{(m_{0}+m_{1})[(ax_{1}+x_{2})^{2}+y_{2}^{2}+(az_{1}+z_{2})^{2}]^{3/2}}\\[8.5359pt] \quad+x_{2}\dot{\theta}^{2}+2\dot{y}_{2}\dot{\theta}+y_{2}\ddot{\theta}\end{array}$$

(7)

$$\begin{array}[]{l}\ddot{y}_{2}=-\frac{mm_{1}y_{2}}{(m_{0}+m_{1})[(x_{1}-x_{2})^{2}+y_{2}^{2}+(z_{1}-z_{2})^{2}]^{3/2}}-\frac{mm_{0}y_{2}}{(m_{0}+m_{1})[(ax_{1}+x_{2})^{2}+y_{2}^{2}+(az_{1}+z_{2})^{2}]^{3/2}}\\[8.5359pt] \quad+y_{2}\dot{\theta}^{2}-2\dot{x}_{2}\dot{\theta}-x_{2}\ddot{\theta}\end{array}$$

(8)

$$\begin{array}[]{l}\ddot{z}_{2}=\frac{mm_{1}(z_{1}-z_{2})}{(m_{0}+m_{1})[(x_{1}-x_{2})^{2}+y_{2}^{2}+(z_{1}-z_{2})^{2}]^{3/2}}-\frac{mm_{0}(ax_{1}+x_{2})}{(m_{0}+m_{1})[(ax_{1}+x_{2})^{2}+y_{2}^{2}+(az_{1}+z_{2})^{2}]^{3/2}}\end{array}$$

(9)

We should remark some essential differences between the restricted and the general 3D-TBP described in the rotating frame. In the general problem, the primaries do not remain fixed on the rotating axis $Gx$, but move on the rotating plane $xz$. Also, the rotating frame $Gxyz$ does not revolve with a constant angular velocity, see Eq. (3), and its origin $G$ does not remain fixed with respect to the inertial system $OXYZ$.

In our numerical computations we should always take that $P_{1}$ is the inner planet, $P_{2}$ is the outer one and the total mass of the system is $m=1$. Also, we should always take that one of the planets (either the inner or the outer) has the mass of Jupiter, namely $m_{i}=10^{-3}$, with either $i=1$ or $i=2$.

Considering the rotating frame and the Poincaré map $y_{2}=0$, periodic orbits are defined as fixed or periodic points of the Poincaré map satisfying the conditions

$$\begin{array}[]{l}x_{1}(0)=x_{1}(T),\;x_{2}(0)=x_{2}(T),\;z_{2}(0)=z_{2}(T),\\ \dot{x}_{1}(0)=\dot{x}_{1}(T),\;\dot{x}_{2}(0)=\dot{x}_{2}(T),\;\dot{z}_{2}(0)=\dot{z}_{2}(T),\dot{y}_{2}(0)=\dot{y}_{2}(T),\end{array}$$

(10)

provided that $y_{2}(0)=y_{2}(T)$ and $T$ is the period. Due to the existence of the Jacobi integral, one of the above conditions is always fulfilled, when the rest ones are fulfilled.

There exist four transformations under which the Lagrangian (1) remains invariant and, subsequently, there exist four symmetries for the orbits of the 3D-GTBP in the rotating frame (Michalodimitrakis 1979a). A periodic orbit defined by (10) is called symmetric, if it is invariant under one of the symmetries of the system. In the planetary TBP, there exist periodic orbits which are symmetric with respect either to the $xz$-plane or to the $x$-axis.

A periodic orbit obeying the symmetry with respect to $xz$-plane has two perpendicular crossings with the $xz$-plane. Planet, $P_{1}$, moves on that plane and, without loss of generality, we may consider that planet, $P_{2}$, is also on the $xz$-plane at $t=0$. Thus, the initial conditions of such a periodic orbit should be

$$\begin{array}[]{llll}x_{1}(0)=x_{10},&\quad x_{2}(0)=x_{20},&\quad y_{2}(0)=0&\quad z_{2}(0)=z_{20},\\ \dot{x}_{1}(0)=0,&\quad\dot{x}_{2}(0)=0,&\quad\dot{y}_{2}(0)=\dot{y}_{20},&\quad\dot{z}_{2}(0)=0.\end{array}$$

(11)

Subsequently, a $xz$-symmetric periodic orbit can be represented by a point in the four-dimensional space of initial conditions

$$\Pi_{4}=\{(x_{10},x_{20},z_{20},\dot{y}_{20})\}.$$

If an orbit has two perpendicular crossings with the $x$-axis, then it is symmetric with respect to that axis. We can assume that at $t=0$ both planets start perpendicularly from the $x$-axis and the periodic orbit has initial conditions

$$\begin{array}[]{llll}x_{1}(0)=x_{10},&\quad x_{2}(0)=x_{20},&\quad y_{2}(0)=0,&\quad z_{2}(0)=0,\\ \dot{x}_{1}(0)=0,&\quad\dot{x}_{2}(0)=0,&\quad\dot{y}_{2}(0)=\dot{y}_{20},&\quad\dot{z}_{2}(0)=\dot{z}_{20}.\end{array}$$

(12)

Thus, a $x$-symmetric periodic orbit can be represented by a point in the four-dimensional space of initial conditions

$$\Pi^{\prime}_{4}=\{(x_{10},x_{20},\dot{y}_{20},\dot{z}_{20})\}.$$

By changing the value of $z_{2}$ (for the $xz$-symmetry) or $\dot{z}_{2}$ (for the $x$-symmetry) a monoparametric family of periodic orbits is formed. Also, a monoparametric family can be formed by changing the mass of a planet, $P_{1}$ or $P_{2}$, but keeping the value $z_{2}$ (or $\dot{z}_{2}$) constant.

The linear stability of the periodic orbits can be found by computing the monodromy matrix $\mathbf{M}$ of the variational equations of the system Eqs. (6)-(9). Since the system is of four degrees of freedom, $\mathbf{M}$ is an $8\times 8$ symplectic matrix and has four pairs of conjugate eigenvalues. One pair is always equal to unity, because of the existence of the energy integral. In order for the periodic orbit to be linearly stable, all of the eigenvalues have to be on the unit circle. Single instability is considered if there exists one pair of real eigenvalues $\{\mu,\mu^{-1}\}$. If an eigenvalue $\mu$ is complex, then there exist also the eigenvalues $\bar{\mu}$, $\mu^{-1}$ and $\bar{\mu}^{-1}$ and the stability is characterized as complex instability, while the remaining eigenvalue pair may either lie or not on the unit circle (Marchal 1990; Hadjidemetriou 2006b). Nevertheless, apart from the pair $(1,1)$, $\mathbf{M}$ may have two or three pairs of real eigenvalues.

When the planetary masses are very small compared to the mass of the Star, periodic orbits in the rotating frame correspond to almost Keplerian planetary orbits in space with orbital elements $a_{i}$ (semimajor axis), $e_{i}$ (eccentricity), $i_{i}$ (inclination), $\omega_{i}$ (argument of pericenter), $\Omega_{i}$ (longitude of ascending node) and $\lambda_{i}$ (mean longitude) or $M_{i}$ (mean anomaly), where $i=1,2$ indicates the particular planet.

A periodic orbit defines an exact resonance for the underlying dynamical system and is locked to a mean motion resonance $n_{P}/n_{J}=(p+q)/p$, where $n_{J}$ is the mean motion of the Jupiter and $n_{P}$ is the mean motion of the other planet. In this paper, we study the $1/2$ (Jupiter is the inner planet $P_{1}$) or the $2/1$ (Jupiter is the outer planet $P_{2}$) resonant cases. We distinguish 2/1 from 1/2 resonance because different dynamics is obtained for the restricted problem. However, this distinction is not essential for the general problem. We can define the resonant angles $\sigma_{1}$ and $\sigma_{2}$ of the system as:

$$\begin{array}[]{c}q\sigma_{1}=(p+q)\lambda_{1}-p\lambda_{2}-q\varpi_{1}\\ q\sigma_{2}=(p+q)\lambda_{1}-p\lambda_{2}-q\varpi_{2}\end{array}$$

In the following, we refer to the resonant angles $\sigma_{1}$, $\Delta\varpi$=$\sigma_{1}-\sigma_{2}$ and $\Delta\Omega=\Omega_{2}-\Omega_{1}$. For the initial conditions (11) the orbital elements $\omega_{i}$ and $\Omega_{i}$ can take the values $\frac{\pi}{2}$ or $\frac{3\pi}{2}$. For the initial conditions (12) the orbital elements $\omega_{i}$ and $\Omega_{i}$ can take the values $0$ or $\pi$. Also, for both symmetries the planets at $t=0$ are at their periastron or apoastron. Subsequently, the resonant angles that correspond to a periodic orbit can take values $0$ or $\pi$.

In Fig. 1, we illustrate an 1/2 resonant stable $xz$-symmetric periodic orbit¹¹1This orbit belongs to the family $F^{1/2}_{g1,i}$, see Section 4.2. in the inertial frame and in the projection space $x_{2}y_{2}z_{2}$ of the rotating frame and for planetary masses $m_{1}=10^{-3}$ and $m_{2}=4\times 10^{-4}$. The initial conditions correspond to the planetary orbital elements

$$\begin{array}[]{cccccc}a_{1}=0.45,&e_{1}=0.66,&i_{1}=12.79^{\circ},&\Omega_{1}=270^{\circ},&\omega_{1}=90^{\circ},&M_{1}=0^{\circ}\\ a_{2}=0.71,&e_{2}=0.34,&i_{2}=20.46^{\circ},&\Omega_{2}=90^{\circ},&\omega_{2}=270^{\circ},&M_{2}=180^{\circ}.\end{array}$$

Its $xz$ symmetry is revealed by the projections of the orbit in the planes of the space $x_{2}y_{2}z_{2}$ (see Fig. 1b). The corresponding resonant angles take the values $\sigma_{1}=0^{\circ}$, $\Delta\varpi=0^{\circ}$ and $\Delta\Omega=180^{\circ}$.

In Fig. 2, we show an 1/2 resonant $x$-symmetric unstable periodic orbit²²2This orbit belongs to the family $G^{1/2}_{g1,i}$, see Section 4.2.. The planetary masses are $m_{1}=10^{-3}$ and $m_{2}=10^{-5}$ and the initial conditions correspond to the planetary orbital elements

$$\begin{array}[]{cccccc}a_{1}=0.33,&e_{1}=0.42,&i_{1}=0.26^{\circ},&\Omega_{1}=0^{\circ},&\omega_{1}=0^{\circ},&M_{1}=0^{\circ}\\ a_{2}=0.52,&e_{2}=0.20,&i_{2}=19.79^{\circ},&\Omega_{2}=0^{\circ},&\omega_{2}=0^{\circ},&M_{2}=180^{\circ}.\end{array}$$

Its $x$ symmetry is revealed by the projections of the orbit in the planes of the space $x_{2}y_{2}z_{2}$ (see Fig. 2b). The corresponding resonant angles take the values $\sigma_{1}=0^{\circ}$, $\Delta\varpi=0^{\circ}$ and $\Delta\Omega=0^{\circ}$.

The distribution of eigenvalues of the linear stability analysis is presented in Fig. 3. Regarding the stable $xz$-symmetric periodic orbit, all eigenvalues lie on the unit circle (Fig. 3a). One pair (of non unit eigenvalues) is very close to 1 and this is shown in the corresponding magnification (Fig. 3c). Concerning the unstable $x$-symmetric orbit there exists one pair of real eigenvalues indicating the instability (Fig. 3b). Also, one pair of eigenvalues is very close to 1 (see the magnification in Fig. 3d).

In Fig.4, we present two different distributions of eigenvalues obtained for two orbits that belong to the family $F^{1/2}_{g2,i}$ (see Section 4.2). In the first case, we have the presence of complex instability, while in the second one, two pairs of real eigenvalues exist. Along a family of periodic orbits the eigenvalues move smoothly either on the unit circle, or the real axis and can either enter the real axis, or the unit circle. Therefore, the stability type can change along the families. In the rest of the paper, we will indicate whether an orbit is stable or unstable without declaring the particular type of instability.

It is known, that in the phase space, a stable periodic orbit is surrounded by invariant tori, where the motion is regular and stays in the neighbourhood of the periodic orbit. In the case of an unstable periodic orbit, chaotic domains exist at least close to the periodic orbit. In general, starting from the neigbourhood of an unstable periodic orbit, the evolution of planetary system shows instabilities and chaotic behaviour. Numerical integrations show that such instabilities appear stronger as far as highly eccentric orbits are concerned (Voyatzis 2008). Hereafter, we present some numerical results that show the above mentioned behaviour of orbits.

We consider the initial conditions of the unstable $x$-symmetric periodic orbit given in the previous section (Fig. 2). Since the initial conditions of the periodic orbit are not exact, the instability leads the evolution far from the periodic orbit. In Fig. 5, we present the evolution of the orbital elements $a_{i}$, $e_{i}$ and $i_{i}$ and the resonant angles $\sigma_{1}$, $\Delta\varpi$ and $\Delta\Omega$. We observe that the eccentricity and the inclination of the heavier planet $P_{1}$ remain almost constant during the whole time of integration. However, the planet $P_{2}$ after the critical time $12\times 10^{3}$ t.u. destabilizes, its eccentricity increases and then oscillates around a high value. Also, its inclination takes large values and shows large oscillation around $90^{\circ}$, i.e. its motion turns from prograde to retrograde and vice versa. A significant oscillation of the semimajor axis of $P_{2}$ after $12\times 10^{3}$ t.u. is also, apparent, though the system remains in the mean motion resonance. The resonant angles $\sigma_{1}$ and $\Delta\varpi$ initially librate arround $0^{\circ}$, while $\Delta\Omega$ increases. After the critical time, $\sigma_{1}$ and $\Delta\varpi$ start to rotate, while $\Delta\Omega$ librates irregularly in the domain $[52^{\circ},145^{\circ}]$ . The computation of the DFLI (see Fig. 8, case II) shows clearly the chaotic nature of the evolution.

If we start with the initial conditions given for the stable $xz$-symmetric periodic orbit (Fig. 1), we will obtain that the orbital elements $a_{i}$, $e_{i}$, $i_{i}$ and the resonant angles remain constant for ever (assuming sufficient integration accuracy). Starting near the periodic orbit, and particularly considering a deviation of $5^{\circ}$ in $\Omega_{1}$ (i.e. we set now $\Omega_{1}=275^{\circ}$), we obtain a regular evolution as it is shown in Fig. 6. The inclination and the eccentricity values of both planets oscillate arround their initial values, while the resonant angles $\sigma_{1}$, $\Delta\varpi$ and $\Delta\Omega$ show librations around the values which correspond to the periodic orbit, namely $0^{\circ}$, $0^{\circ}$ and $180^{\circ}$, respectively. Now, the DFLI does not increase (Fig. 8, case I)

If we start not sufficiently close to the stable periodic orbit, the stability of evolution is not guaranteed. We consider a larger deviation, compared to the deviation used in the previous case. Particularly, we consider a deviation of $15^{\circ}$ in $\Omega_{1}$ and we integrate the orbit. The results are shown in Fig. 7. The chaotic behaviour of the evolution is obvious and finally the planet $P_{2}$ is scattered. The resonant angles $\sigma_{1}$ and $\Delta\varpi$ rotate, but it is remarkable that $\Delta\Omega$ librates irregularly around $180^{\circ}$. The DFLI clearly shows the chaotic evolution (Fig. 8, case III).

We consider the restricted problem with Jupiter moving in a circular planar motion with period $2\pi$ and a massless planet which moves in the space. In this problem, we can obtain families of periodic orbits which bifurcate from vertical critical periodic orbits of the planar restricted problem (Hénon 1973). These families exist either when Jupiter is the inner, or the outer planet (resonance either 1/2, or 2/1, respectively).

Following the notation of Kotoulas and Voyatzis (2005), families of periodic orbits which are symmetric with respect to the $xz$-plane are denoted by $F$ and families of periodic orbits symmetric with respect to the $x$-axis are denoted by $G$. The subscript $c$ or $e$ indicates the spatial restricted problem (circular or elliptic, respectively) from which these families bifurcate. Moreover, $m$ indicates the continuation with respect to the mass of the initially massless body. Also, the corresponding resonance is indicated by a superscript.

We computed and present the families of 2/1 and 1/2 resonant symmetric periodic orbits of the 3D-RTBP in Fig. 9 (see also Hadjidemetriou and Voyatzis (2000) and Kotoulas (2005), respectively). We denote these families by $F_{c,i}^{2/1}$, where the subscript $i$ indicates the continuation from plane to space (see Section 4.2). For the 2/1 resonance such families correspond to relatively high eccentricity values of the massless planet $P_{1}$. Conversely, in the 1/2 resonance all periodic orbits correspond to almost circular orbits of the massless planet $P_{2}$ and are unstable.

It is known that the periodic orbits of the 3D-RCTBP of period $T$ can be analytically continued, by increasing the planetary mass of the initially massless body to a periodic orbit of the 3D-GTBP of the same symmetry, provided that $T\neq 2k\pi$, where $k$ is an integer (Ichtiaroglou and Michalodimitrakis 1980). This is the case for all orbits of the families given in Fig. 9. For the particular computations presented hereafter, we continue analytically the periodic orbits indicated by the dots in Fig. 9 with respect to the mass $m_{1}$ (inner planet) or $m_{2}$ (outer planet) according to the resonance 2/1 or 1/2, respectively. These periodic orbits correspond to some initial inclination value of the massless body, $i_{10}$ or $i_{20}$, and identify the generated families.

Resonant symmetric periodic orbits for the planar general problem (2D-GTBP) can be obtained by starting from either the circular family of orbits (e.g. Hadjidemetriou (2006a)), or from resonant families of the planar restricted problem (e.g. Antoniadou et al. (2011)). Particularly, for the 2/1 (or 1/2) resonance, planar periodic orbits have been computed by (Beaugé et al. 2006; Voyatzis and Hadjidemetriou 2005; Voyatzis et al. 2009). Now, the distinction between the 2/1 and 1/2 resonance is not essential, since both planetary masses are nonzero and we refer to the mass ratio $\rho=m_{2}/m_{1}$, where $P_{1}$ is the inner planet. In the following computations, we always set $m_{1}=0.001$ (Jupiter’s mass). In Fig. 16, we present the 1/2 resonant families $f_{1}$ and $f_{2}$ of the general planar problem. The families are given for various values of the mass ratio $\rho$ in the projection space defined by the planetary eccentricities.

The variational equation that corresponds to the $z$-component of the motion of planet $P_{2}$ is written as³³3Note that the variational equation given in (Michalodimitrakis 1979a) holds only for the normalization $m_{0}+m_{1}=1$.

$$\ddot{\zeta}_{2}=A\zeta_{2}+B\dot{\zeta}_{2}$$

(13)

where

$$\begin{array}[]{l}A=-\frac{mm_{0}[1-\frac{m_{2}(\dot{\theta}x_{2}+\dot{y}_{2})}{m\dot{\theta}x_{1}}]}{(m_{0}+m_{1})[(\frac{m_{1}x_{1}}{m_{0}}+x_{2})^{2}+y_{2}^{2}]^{3/2}}-\frac{mm_{1}[1+\frac{m_{0}m_{2}(\dot{\theta}x_{2}+\dot{y}_{2})}{mm_{1}\dot{\theta}x_{1}}]}{(m_{0}+m_{1})[(x_{1}-x_{2})^{2}+y_{2}^{2}]^{3/2}}\\[8.5359pt] B=\frac{m_{0}m_{2}y_{2}}{(m_{0}+m_{1})\dot{\theta}x_{1}[(x_{1}-x_{2})^{2}+y_{2}^{2}]^{3/2}}-\frac{m_{0}m_{2}y_{2}}{(m_{0}+m_{1})\dot{\theta}x_{1}[(\frac{m_{1}x_{1}}{m_{0}}+x_{2})^{2}+y_{2}^{2}]^{3/2}}.\end{array}$$

(14)

If $\Delta(T)=\{\xi_{ij}\}$, $i,j=1,2$, is the monodromy matrix of Eq. 13 for a planar periodic orbit of period $T$, then, this orbit has a vertical critical index (Michalodimitrakis 1979a) equal to

$$a_{v}=\frac{\xi_{11}\xi_{22}+\xi_{21}\xi_{12}}{\xi_{11}\xi_{22}-\xi_{21}\xi_{12}},$$

(15)

Any periodic orbit of the planar problem with $|a_{v}|=1$ is vertical critical (v.c.o.) and can be used as a starting orbit for analytical continuation in the spatial problem with the inclination of $P_{2}$ as a parameter (in computations we increase $z_{2}$ or $\dot{z}_{2}$).

In Fig. 16, we present the v.c.o. by dots (see also Tables 1-3). In Fig. 17a, we illustrate the variation of the index $a_{v}$ along the family $f_{1}$ for some values of $\rho$. For $\rho<0.0372$, we find three v.c.o. and starting from one of them (see Table 1) we can analytically continue $xz$-symmetric periodic orbits in space. We denote such families by $F^{1/2}_{g1,i}$ where the subscript $g$ indicates that the family bifurcates from the general planar problem (it is followed by a number that states the generating planar family) and $i$ indicates that the continuation takes place from plane to space. From the remaining couple of v.c.o. (see Table 2) $x$-symmetric periodic orbits in space are obtained forming the families $G^{1/2}_{g1,i}$. An example of the above families is given in Fig. 18a in the projection space $e_{1}e_{2}\Delta i$, where $\Delta i$ is the mutual inclination. We obtain that the family $G^{1/2}_{g1,i}$ starts from and ends at v.c.o. of the same planar family $f_{1}$. For $\rho>0.0372$ the couple of v.c.o., from which we obtained the family $G^{1/2}_{g1,i}$, disappears. This transition is explained by the variation of the index $a_{v}$ along the family $f_{1}$ and is shown in Fig. 17b.

Concerning the families $F^{1/2}_{g1,i}$, up to the mass ratio value $\rho^{*}\approx 0.12$ all the corresponding v.c.o. are unstable and, subsequently, the continuation generates families that start as unstable. The families continue to be unstable until the end of the computations. For $\rho>\rho^{*}$ the v.c.o. belong to the stable segment of family $f_{1}$. Now the generated families $F^{1/2}_{g1,i}$ start as stable and show a monotonical increase of the inclinations $i_{1}$ and $i_{2}$, and, subsequently, $\Delta i$, too (see Fig. 18). The resonant angles that correspond to these periodic orbits are

$$\Delta\varpi=0^{\circ},\;\Delta\Omega=180^{\circ},\;\sigma_{1}=0^{\circ}.$$

The stability type along the families $F^{1/2}_{g1,i}$ changes above a maximum value for the mutual inclination, $\Delta i_{max}$ that depends on $\rho$, as it is shown in Fig. 19. In case (a), Jupiter is the inner planet, $P_{1}$, and $P_{2}$ takes large mass values (up to 20 Jupiter masses), as $\rho$ increases. In case (b), we present the same computations, but now Jupiter is the outer planet, $P_{2}$ and the other planet, $P_{1}$, takes smaller and smaller mass’ values, as $\rho$ increases. We observe that we can have stable orbits up to mutual inclinations $40^{\circ}<\Delta i_{max}<50^{\circ}$ and there is no essential difference, whether Jupiter is the inner, or the outer planet in the system.

In Fig. 20, we present the families $F^{1/2}_{g2,i}$, which bifurcate from the v.c.o. of families $f_{2}$ given in Fig. 16b. The corresponding resonant angles are

$$\Delta\varpi=180^{\circ},\;\Delta\Omega=180^{\circ},\;\sigma_{1}=180^{\circ}.$$

The families which start from stable v.c.o. (in the range $0.3<\rho<7$, approximately) show a segment of stable periodic orbits up to a critical mutual inclination value, $\Delta i_{max}$. $\Delta i_{max}$ seems to increase, as $\rho$ increases in the above mentioned interval and reaches $30^{\circ}$, approximately. The unstable segments show transitions from real to complex instability and vice versa (see e.g. Fig. 4).