# Sitnikov problem

Post-publication activity

Curator: Rudolf Dvorak

The Sitnikov problem is a special case of the restricted three-body problem that allows oscillatory motions (Sitnikov 1961).  Figure 1: animation of the motions in the MacMillan problem. The parameters are $e=0$, the initial conditions are $z(0)=1$ and $\dot z(0)=0$.  Figure 2: animation of the motions in the Sitnikov problem. The parameters are $e=0.1$, the initial conditions are $z(0)=1$ and $\dot z(0)=0$.

## Introduction

The Sitnikov problem is a sub-case of the spatial elliptic restricted three-body problem that allows oscillatory type of motions: a massless body moves (oscillates) along a straight line that is perpendicular to the orbital plane that is formed by two equally massed primary bodies moving on symmetric Keplerian orbits (cf. Figure 1,2). In its original formulation (Sitnikov 1961) the equation of motion for the massless body turns out to be $1 \frac{1}{2}$ dimensional and depends only on time $t$ and on one parameter, the eccentricity of the primaries $e$. The problem is a direct generalization of the so-called MacMillan problem (MacMillan 1911) that is the integrable approximation of the Sitnikov problem for $e=0$. For $e\neq0$ the dynamical system is chaotic, despite its simple form, and allows all different kinds of motions that are found in chaotic dynamical systems. The simplicity of its mathematical formulation together with the complexity of its possible kinds of motions make the Sitnikov problem a unique dynamical problem in the context of celestial mechanics.

## Derivation from Newton's law of gravity

Let $m_1$, $m_2$, $m_3$ be the masses of three spherical particles $P_1$, $P_2$, $P_3$, where $P_1$, $P_2$ are called primaries, and $P_3$ is called the third body. We assume an inertial coordinate system with the origin in the center of gravity. Let us denote in this reference system by $\vec r_i$ the position vector of the particle $P_i$, and by $\Delta_{i,j}$ the scalar distance between the particles $P_i$ and $P_j$, with $i\neq j$ and $i,j=1,2,3$. The gravitational force $\vec F_{ij}$ acting on one of the particles $P_i$ of mass $m_i$ due to the particle $P_j$ of mass $m_j$ is then given by:

\begin{equation} \vec F_{ij}=-G \sum _{i\neq j,j=1}^3 m_i m_j\frac{ \left(\vec r_i-\vec r_j\right)}{\Delta_{i,j}^3} \ , \end{equation}

where $G$ is the gravitational constant. Let $\vec r_k=(x_k,y_k,z_k)$ with $k=1,2,3$. The Euclidean distance $\Delta_{i,j}$ is then given by:

\begin{equation} \Delta _{i,j}=\sqrt{ \left(x_i-x_j\right)^2 +\left(y_i-y_j\right)^2 +\left(z_i-z_j\right)^2} \ , \end{equation}

If we denote the accelerations by $\ddot{\vec r_k}=(\ddot x_k, \ddot y_k, \ddot z_k)$ with $k=1,2,3$, then the system of equations of motions is given by:

\begin{equation} \ddot{\vec r_1}=-G \left(\frac{m_2 \left(\vec r_1-\vec r_2\right)}{\Delta _{1,2}^3} +\frac{m_3 \left(\vec r_1-\vec r_3\right)}{\Delta _{1,3}^3}\right) \end{equation}

\begin{equation} \ddot{\vec r_2}=-G \left(\frac{m_1 \left(\vec r_2-\vec r_1\right)}{\Delta _{1,2}^3} +\frac{m_3 \left(\vec r_2-\vec r_3\right)}{\Delta _{2,3}^3}\right) \end{equation}

\begin{equation} \ddot{\vec r_3}=-G \left(\frac{m_1 \left(\vec r_3-r_1\right)}{\Delta _{1,3}^3} +\frac{m_2 \left(\vec r_3-\vec r_2\right)}{\Delta _{2,3}^3}\right) \ . \end{equation}

The above system can be solved together with initial conditions that are given at initial time $t=t_0$:

\begin{equation} \vec r_k(t_0)=(x_k(t_0), y_k(t_0), z_k(t_0)) \ , \ \dot{\vec r_k}(t_0)=(\dot x_k(t_0), \dot y_k(t_0), \dot z_k(t_0)) \ \end{equation}

with $k=1,2,3$.

We obtain the equation of motion of the classical Sitnikov problem from the above system by choosing a proper orientation of the basic reference axes, by making additional assumptions on the involved masses, and by restricting the motions of $P_1$, $P_2$, $P_3$ in terms of special initial conditions. As we will show in a moment, with this it will become possible to reduce the number of degrees of freedom of the above system of 9 autonomous differential equations of second order to 1 (!) single non-autonomous equation of motion, that defines the so called Sitnikov problem. Let $(X,Y,Z)$ denote the axes of our inertial coordinate system.

i) We choose the basic plane of reference ($X$, $Y$) to coincide with the orbital planes of $P_1$, $P_2$. Let $m_1=m_2\equiv m$ be of the same mass. The motion of the primaries, in this plane, in absence of the mass $m_3$ takes place along two antisymmetric Keplerian ellipses with a common focus that coincides with the center of gravity of this sub-system. Let, in addition, the $X$-axis be parallel to the line that connects the two apocenters of these ellipses. Then, the anti-symmetry between the unperturbed motions of $P_1$, $P_2$ is such that

\begin{equation} \vec r_1=-\vec r_2=(x_1,y_1,0) \ , \end{equation}

where we have already set $z_1=0$ that is consistent with the assumptions that we make on the third body:

ii) Let us assume that 1) the mass of the third body can be neglected with respect to the masses of the primary bodies, and 2) that the motion of the third body starts along the $Z$ axis with special initial conditions, such that the velocity components of $P_3$ that are normal to the $Z$-axis vanish:

We formally set $m_3=0$, $x_3(0)=y_3(0)=0$, $\dot x_3(0)=\dot y_3(0)=0$ and thus find $\ddot x_3=\ddot y_3=0$.

The assumptions we make allow us to investigate the reduced set of differential equations:

\begin{equation} \ddot{x_1}=-\frac{G m x_1}{4\left(x_1^2+y_1^2\right)^{3/2}} \end{equation}

\begin{equation} \ddot{y_1}=-\frac{G m y_1}{4\left(x_1^2+y_1^2\right)^{3/2}} \end{equation}

\begin{equation} \ddot{z_3}=-\frac{2Gm z_3}{\left(x_1^2+y_1^2+z_3^2\right)^{3/2}} \ . \end{equation}

We notice that the first two of the above equations are uncoupled from the third. Moreover, the coupling terms in the third equation are of the special form

\begin{equation} r_1^2=x_1^2+y_1^2 \ , \end{equation}

that is the distance of one of the primaries from the common barycenter of the system. Let us drop the index $1$ in $r_1$ from now on. The equations of motion of the unperturbed Kepler problem in polar coordinates are given by:

\begin{equation} \ddot r-r\dot\phi^2=-\frac{\mu}{r^2} \ , \ \end{equation}

\begin{equation} \frac{1}{r}\frac{d}{dt}\bigg(r^2\dot\phi\bigg)=0 , \end{equation}

where

\begin{equation} \mu=G(m_1+m_2)=2Gm \end{equation}

is the mass parameter, and the angular variable $\phi$ is called the true anomaly. We make use of Kepler's second law, and the fact that the angular momentum $c$ is a constant of motion that depends on the masses, the semi-major axes $a$, and the eccentricities of the primary bodies $e$, see Dvorak & Lhotka (2013):

\begin{equation} c=r^2\dot\phi=\sqrt{\mu a(1-e^2)} \ . \end{equation}

We are thus able to express the radial component of the unperturbed two-body motion in terms of $\mu$, $c$ only:

\begin{equation} \ddot r-r\bigg(\frac{c}{r^2}\bigg)^2=-\frac{\mu}{r^2} \ . \end{equation}

We drop the index $3$ in $z_3$ in the remaining equation of motion for the third body, and we are left with the equations:

\begin{equation} \ddot z=-\frac{2Gmz}{\big(r^2+z^2\big)^{3/2}} \ , \ \end{equation}

\begin{equation} \ddot r = -\frac{2Gm}{r^2}+\frac{c^2}{r^3} \ . \end{equation}

The solution for $r$ is given by:

\begin{equation} r(t)=a(1-e \cos (E(t))) \end{equation}

where $E=E(t)$ is the eccentric anomaly of one of the primary bodies, that is related to the mean motion $n$ and mean anomaly $M=M(t)=n(t-t_0)$ through Kepler's equation: \begin{equation} M(t)=E(t)-e \sin (E(t)) \ . \end{equation}

We notice that by Kepler's second law we have $n^2a^3=\mu$. If we then substitute for $r=r(t)$ we are left with:

\begin{equation} \ddot{z}+\frac{2Gmz}{\left(r(t)^2+z^2\right)^{3/2}}=0 \ , \end{equation}

that is a $1 \frac{1}{2}$ dimensional, non-autonomous equation of motion that describes the restricted motion of the third body along the $Z$-axis under our assumptions that we made.

If we denote by $T$ and $V$ the kinetic and potential energies of the system, then the dynamical equation may also be derived from a Lagrangian

\begin{equation} L=T-V=\frac{\dot z^2}{2}+\frac{2Gm}{\sqrt{r(t)^2+z^2}} \ , \end{equation}

using

\begin{equation} \frac{d}{dt}\bigg(\frac{\partial L}{\partial\dot z}\bigg)-\frac{\partial L}{\partial z}=0 \end{equation}

or from a Hamiltonian:

\begin{equation} H=T+V=\frac{v^2}{2}-\frac{2Gm}{\sqrt{r(t)^2+z^2}} \ , \end{equation}

using

\begin{equation} \frac{dz}{dt}=\frac{\partial H}{\partial v} = v = \dot z \end{equation}

\begin{equation} \frac{dv}{dt}=-\frac{\partial H}{\partial z} = - \frac{2Gmz}{\left(r(t)^2+z^2\right)^{3/2}} \end{equation}

(where we used $v=\dot z$). We notice that by a proper choice of the parameters and units we are able to simplify the problem further. We set the total mass of the system $2Gm=G(m_1+m_2)=1$, and set the unit of distance $a=1$ (we notice that for some studies it is convenient to set the distance between the two primaries $2a=1$). Then, from Kepler's second law $n^2a^3=1$, we find the unit of time in terms of the mean motion $n=1$ such that the period, $T=2\pi/n$, of one revolution of the primaries takes the value $T=2\pi$. Moreover, if we set $t_0=0$, then we find $M=t$ and that the primaries start at their pericenter at $t=0$.

## The MacMillan problem

There exists an integrable case of the Sitnikov problem for $e=0$. This problem is usually called the MacMillan problem (Pavanini 1907, MacMillan 1911). In normalized units it is defined by the Hamiltonian

\begin{equation} H=\frac{v^2}{2}-\frac{1}{\sqrt{1+z^2}} \end{equation}

that results in the time independent equation of motion

\begin{equation} \ddot z + \frac{z}{\big(1+z^2\big)^{3/2}}=0 \ . \end{equation}

The direct integration with respect to time gives:

\begin{equation} \bigg(\frac{dz}{dt}\bigg)^2=\frac{2}{\sqrt{1+z^2}}-2C \ , \end{equation}

where $C$ is a constant that depends on the initial conditions. It is shown in MacMillan 1911 that by the substitutions

\begin{equation} 1+z^2=\frac{1}{u^2} \end{equation}

and

\begin{equation} v^2=\frac{1-u}{1-C} \end{equation}

the resulting integral

\begin{equation} \int_0^tdt=\int_0^v\frac{dv}{(1-2k^2v^2)^2\sqrt{(1-v^2)(1-k^2v^2)}} \end{equation}

can be solved in terms of elliptic functions, where

\begin{equation} k^2=\frac{1}{2}(1-C) \ . \end{equation}

The motion of the third body along the $Z$-axis is periodic provided that $0<C\leq1$. For small $k^2<1/2$ the period $P$ in power series of $k^2$ is then given by:

\begin{equation} P=2\pi\bigg(1+\frac{9}{4}k^2+\frac{345}{64}k^4+\frac{3185}{256}k^6+ \frac{457065}{16384}k^8+\dots\bigg) \end{equation}

The integrability of the system for $e=0$ has been shown in Pavanini (1907), MacMillan (1911). A modern treatment of the problem in terms of elliptic functions can be found in Belbruno et al. (1994).

## Oscillatory Motion in the Restricted Three-Body Problem

There exists bounded and unbounded motions in the Sitnikov and MacMillan problem. If $|z(0)|+|\dot z(0)|$ is small enough then $z(t)$ remains bounded for all time. If $|\dot z(0)|$ is large enough then $z(t)\to\infty$ for $t\to\infty$. In the case $e\neq0$ the motion close to the escape velocity, say $\dot z_{crit}$, is of complex nature due to the non-integrability of the dynamical system.

In Sitnikov (1961) the author showed for the first time the existence of oscillatory motions in the restricted three-body problem that are unbounded. This exciting result is a consequence of the following theorem (after Moser 1973):

Theorem Given a sufficiently small eccentricity $e>0$ there exists an integer $m=m(e)$ such that any sequence $s$ with $s_k\geq m$ corresponds to a solution of the differential equation of the Sitnikov problem.

Here, $s=[\dots,s_{k-1},s_k,s_{k+1},\dots]$ is a doubly infinity sequence of integers $s_k$ with $s_k\geq m$, and the $s_k$ are defined by

\begin{equation} s_k=\bigg[\frac{t_{k+1}-t_k}{2\pi}\bigg] \ , \end{equation}

where $[.]$ denotes the integer part, and $t_k$ are discrete times at which $z(t_k)=0$, with $t_k<t_{k+1}$. The $s_k$ are therefore a measure of the number of complete revolutions of the primaries between two consecutive zeroes of the solution $z(t)$ (see Figs. 5-7). We remark that the $s_k$ can be chosen completely independently! Thus, if one chooses an unbounded sequence of integers $s_k$ the corresponding solution will be unbounded too, but will contain an infinite number of zeroes.

The result has originally been shown for zero mass of the third body $m_3$, and small $e$. It has been generalized to non-zero mass of the third body $m_3\neq0$, and for the whole eccentricity interval $0<e<1$ in Alekseev (1968a, 1968b, 1969). A concise review of the techniques that are needed to prove these kinds of theorems can be found in Moser (1973). The basic idea is to show the topological equivalence of a discrete mapping, that describes the motion of the Sitnikov problem close to the critical velocity boundary, to the Bernoulli shift, that can itself be studied using symbolic dynamics and the shift operator.

We show some typical kinds of orbits of the Sitnikov problem in Figs.5-7. For $e=0$ (the MacMillan problem) we find an integrable orbit with period $3$ (in terms of revolution periods of the primaries). The initial conditions are $z(0)=1$, $\dot z(0)=0$. In the middle plot we show an orbit for $e=0.3$, $z(0)=2.0$, $\dot z(0)=0$ that corresponds to the integer sequence $7-6-5-\dots$; in the lower plot we show an orbit with $e=0.3$, $z(0)=0.1$, $\dot z(0)=2.5$ that has the starting sequence $8-21-4-4-19-\dots$. We notice that the latter plot shows a chaotic orbit.

## Vrabec-map

To visualize the theorem of Moser, namely whatever sequence of integer numbers $s_k$, which corresponds to revolutions of the primary bodies, one can imagine, there is always one 'real' orbit having these successive crossings of the barycentre (with the restriction $s_k \geq m(e)$).

Following an original idea of F. Vrabec (Dvorak et al. 1993) one can visualize the complexity of the parameter space of the Sitnikov problem in the following way: for a fixed value of the eccentricity $e$ one integrates orbits of the third body -- for a certain number of crossings of the barycenter -- using a grid of initial conditions in true anomaly $0^{\circ} < \phi < 360^{\circ}$ of the primaries versus the velocity $\dot z > 0$ with ($z = 0$) as starting point. From one crossing of the massless body through the barycenter we count the number of full orbits of the primary bodies which gives the number of the Moser sequence. One can characterize these numbers with the aid of a color code: after the first crossing we proceed to the next crossing, count the number of revolutions between the last crossing and this one which now defines the numbers in the Moser sequence shown as the regions in different colors in the diagram. In Fig. 8 one can see a fractal structure of the orbits in the Sitnikov problem after the 6th crossing (after Lang 2011). Figure 8: parameter space study (initial true anomaly vs. velocity) of stable and unstable motion in the Sitnikov-problem. The color code counts the number of crossings per revolution period of the primaries.

## Phase space structure

We show a typical phase portrait in the $(z,\dot z)$-plane for $e=0.1$ in Fig. 9. The central equilibrium is situated at $z=\dot z=0$ and corresponds to the kind of motion where the massless body is at rest. The center is surrounded by KAM curves of motion that form elongated ellipses with their major axis situated along the $\dot z$-axis. At larger distances from the barycenter the closed invariant curves open up, and are replaced by librational islands that are surrounded by chaotic motion that is bounded by further closed invariant curves. The outermost invariant curve (cantorus) has the shape of a diamond and separates the 2:1 resonant island from the central dynamics: it is situated in a chaotic sea that is bounded by the last invariant curve that is elongated to $\pm\infty$. Starting in the white regime outside the outermost invariant curve motion is unbound - the massless particle leaves the system forever. With increasing $e$ more and more invariant curves get destroyed, and higher order islands show up. We show the evolution of the 2:1 resonant island in the animated pictures in Fig. 10 and 11. In Fig. 10 we see the shrinking of the closed invariant curves with increasing $e$, and the presence of more and more chaos, until no stable motion is left (around $e=0.25$). In Fig. 11 we see that the central island of the 2:1 resonant island can reappear for larger values of $e$. However, higher order islands, and the presence of chaos destabilizes the motion close to the resonant island further.

## Linear stability, and approximate nonlinear solution Figure 12: transfer matrix $\rm{Tr}(R)$ versus $e$ obtained from the linearized equations of motion close to the barycenter.

The analysis of the eigenvalues of the linearized equations of motion close to the central equilibrium shows that it is a center. However, the detailed investigation of the trace of the transfer matrix that can be obtained using Floquet theory (Hagel & Lhotka 2005) reveals that the trace of the monodromy matrix may reach the critical values $\pm2$, at which slight perturbations may already harm the stability of the motion. The trace of the transfer matrix, up to order $2$ in eccentricity, is given by (Lhotka 2004):

\begin{equation} \rm{Tr}(R)= 2 \cos \left(\frac{21 e^2}{31}+4\right) \ . \end{equation}

We plot $Tr(R)$ versus $e$ (negative values correspond to the case, where the primaries start at their apocenter) in Fig. 12. Critical values are found close to $\pm0.544\dots$, $\pm0.859$, $\pm0.966$. In Martinez & Chiralt (1993) they authors find two sequences of critical values of $e$, and in addition values at which the central equilibrium becomes unstable ($e=0.8558625$). We numerically integrate such a case in Fig. 13. Figure 13: at $e=0.8558625$ the central equilibrium bifurcates, and the barycenter becomes unstable (after Dvorak, Lhotka 2013).

Advanced perturbative methods may be used to obtain an approximate nonlinear solution of the problem. Using Floquet theory, Courant Snyder transformation, and Poincaré-Lindstedt methods the first order nonlinear solution, for $\dot z=0$ takes the form (Lhotka 2004):

\begin{equation} z(t)=C_1 \cos \left(2 \sqrt{2} \sigma(t) \right) \ , \quad \sigma(t)= t + \frac{48}{31} e \sin (t) \ , \quad C_1=\left(1+\frac{24 e}{31}\right)z(0) \ . \end{equation}

The solution may be used to reproduce the regular kind of motions for small $e$, and small $z(0)$. Higher order expressions of the transfer matrix, and the nonlinear solution may be found in Hagel & Lhotka (2005), Lhotka (2004). A $7$th order solution reproduces the orbit quite well, we show an example orbit for $e=0.1$, $z(0)=0.01$, $\dot z(0)=0$ in Fig. 14,15. At the beginning the two curves (dashed, blue=numerical, black, thick=analytical) overlap completely (Fig. 14), after a long enough time the approximation error in the phase shift accumulates (Fig. 15). Figure 14: comparison of a 7th order analytical and numerical solution for $e=0.1$, $z(0)=0.01$, $\dot z(0)=0$. Figure 15: comparison of a 7th order analytical and numerical solution for $e=0.1$, $z(0)=0.01$, $\dot z(0)=0$.

## Parameter study: the Crab-diagram

The phase space from e=0 to e=0.3

Fig. 16 serves to describe the motion of the planet for values of $0.5 < z_{ini} < 3$ and different eccentricities $0.0 < e < 0.3$ and can be regarded as a synopsis of the former shown surfaces of section. The yellow region characterizes stable motion, the blue region leads to an escape in relatively short time and the red regions on the edge of the yellow one are the so-called sticky orbits (Dvorak et al. 1998) which escape after a very long time of seemingly stable motion. The 'main land' of stable orbits around the barycenter is diminishing from $z_{ini}=1.5$ on to $z_{ini} \sim 1$ with increasing $e$. For small eccentricities the yellow region between $1.5 < z_{ini} < 2.5$ is close to a special periodic orbit which belongs to the 1:2 mean motion resonance, where the primaries fulfill two complete periods whereas the planet makes just one oscillation. One can call the region the '1:2 island'; higher order periodic orbits due to higher resonance in the sense described above appear in form of spikes on the edge of the yellow region with increasing eccentricity. Whenever such an 'island' separates from the big 1:2 island, it itself suffers from a sudden decrease and finally turns into an unstable orbit (blue region). Close to $e \simeq 0.225$ the 1:2 island disappears almost completely, but then reappears again and its size increases. Figure 16: parameter study (eccentricity vs. initial distance); yellow stands for stable orbits, blue for unstable ones.

## Motion of the planet off the z-axis

Using the full set of equations of motion for the third body it is possible to investigate the different kinds of motions off the z-axis also.

In the circular problem an infinite sequence of stability intervals on the z-axis with growing width for larger distances from the barycenter was found using an analytical approach. The first interval on the z-axis is between $5.044 < z < 5.052$; from Table 1 it is evident that the width decreases with the initial distance $z$ and seems to converge to a maximum value.

Table 1: stable intervals on the z-axis
interval low. lim. upp. lim. width
1 5.044 5.052 0.008
2 5.453 5.470 0.017
3 5.848 5.872 0.024
4 6.230 6.259 0.029
5 6.600 6.634 0.034
6 6.961 6.999 0.037
7 7.314 7.353 0.039
8 7.657 7.699 0.042
.. .. .. ..
12 8.965 9.010 0.045
.. .. .. ..
24 12.425 12.475 0.050
25 12.691 12.741 0.050

This is very interesting because even off the z-axis the planet's motion is stable! There exist up to now no analytical approach to determine the radius of such small 'allowed' deviations so that the motion is stable. The numerical simulations can answer this question only in part, but it is evident that these regions are very small compared to the distance from the barycenter. We show an orbit that started inside the stable interval number 12 (see Tab. 1) in Fig. 17. This double bouquet orbit moves inside a double cone extending symmetrically from the $x-y$ plane for $z=0$. Note the different scales for the plane x-y which is 5 orders of magnitude smaller than in the z-direction. An interesting study would be to determine the volume of the double-cone where inside stable orbits exist: is it decreasing or increasing with larger with $z_{\mathrm{ini}}$?

## Literature review

The Sitnikov problem has been generalized or extended in various ways. The system has been studied for unequal masses of the primaries $m_1\neq m_2$, (Perdios & Markelos 1987, Perdios & Markelos 1988), and $m_3\neq0$ (Dvorak & Sun 1997), the spherical point masses have been replaced by oblate bodies also (Pandey & Ahmad 2013, Douskos et al. 2012, Kalantonis et al. 2008). The motion of the primaries out of the $(X,Y)$-plane (Dvorak & Sun 1997, Perdios & Markelos 1988, Perdios & Markelos 1987), the motion of the third body off the $z$-axis has been studied as well (Perdios & Kalantonis 2012, Bountis & Papadakis 2009, Soulis et al. 2008, Soulis et al. 2007, Belbruno et al. 1994), the number of primaries has been increased (Bountis & Papadakis 2009, Soulis et al. 2008), to a continuous ring in Bountis & Papadakis (2009). Additionally, also relativistic forces have been included (Kovács et al. 2011), a complex 'zoo' of very fundamental dynamical properties have been found and investigated: sticky and chaotic orbits (Kovács & Érdi 2009, Hevia & Ranada 1996), escape and diffusion channels (Dvorak 2007), the fractal structure of the phase and parameter space (Kovács & Érdi 2007, Dvorak 1993), bifurcations & families of periodic orbits (Perdios 2007, Perdios & Kalantonis 2006) as well. Various advanced numerical, and analytical methods have been developed and tested on this simple dynamical system: chaos indicators, Poincaré surface of sections and stroboscopic maps (Kovács & Érdi 2007, Dvorak 2007, Jiménez & Escalona-Buendía 2001, Liu et al. 1991), symplectic mappings (Liu & Sun 1990), specialized numerical integration techniques and shadowing orbits (Urminsky 2010), perturbation theory based on normal form theory (Di Ruzza & Lhotka 2011), Floquet theory and Courant Snyder theory (Hagel 2009, Hagel & Lhotka 2005, Faruque 2003, Faruque 2002, Jalali & Pourtakdoust 1997, Wodnar 1995, Hagel Trenkler 1993, Hagel 1992), iterative methods, to name a few. Mathematical theories have been used to prove the existence of special kinds of orbits (Corbera & Llibre 2000, Chesley 1999A, Martinez & Chiralt 1993), linear and non-linear stability theorems have been found (Kalas & Krasil'Nikov 2011, Sidorenko 2011) for different regimes of the phase and parameter space.

Outline of the research studies (reversed chronological order):

In Pandey & Ahmad (2013) the authors study the motions of the third body in the case where the primaries are three oblate spheroids, situated at the edges of an equilateral triangle, in the circular problem. They establish a relation between the oblateness parameter of the primaries and the length of the sides of the equilateral triangle, and also provide the stability region as well as periodic orbits in the planar and 3D problem. []

Families of 3-dimensional periodic orbits, that bifurcate from self-resonant orbits of the Sitnikov family, are studied in Perdios & Kalantonis (2012). It is shown there that the islands of stability of these families disappear gradually in 3-dimensions by varying the mass parameter of the system. []

In a study of the restricted three-body problem with equal prolate primaries the authors of Douskos et al. (2012) show the existence of new equilibrium points that are located along the z-axis and give rise to new type of straight-line periodic oscillations, that are different from the classical motions in the Sitnikov problem. []

The known result, that in a first order approximation the Sitnikov problem has the form of a Hill-type equation (linear second order equation with time dependent, periodic coefficients), is used in Kalas & Krasil'Nikov (2011) to show that the stability of the center equilibrium solution depends on the eccentricity of the primaries $e$. The authors find that the center is stable for almost any $e$, with the exception of a discrete set of $e$ values that accumulates at $e=1$. []

The authors of Di Ruzza & Lhotka (2011) implement a high order (40) Birkhoff normal form to show the nonlinear stability of a small region around the central equilibrium point for long times. The result is based on analytical estimates of the remainder of the normal form of a simplified model that is obtained using perturbation theory. []

Relativistic effects are investigated in Kovács et al. (2011). The study includes the description of the phase-space geometry of the Sitnikov problem in the first post-Newtonian approximation. It is shown that the topology of the phase space strongly depends on the strength of the relativistic pericentre advance. Transient chaos and bifurcations are found, the hyperbolicity of a chaotic saddle is shown numerically.[]

Sidorenko (2011) studies the stability of the vertical periodic motions of the third body in the linear approximation of the classical problem. The author provides a detailed study of the alternation of the stability and instability, within the family of periodic orbits, in dependency of the amplitude of vertical oscillations. []

Shadowing orbits provide a mechanism to show that the dynamics that is represented by numerical simulations that include errors can still be representative for the true dynamics. In Urminsky (2010) the author demonstrates on the basis of the Sitnikov problem, that unstable orbits of the three-body problem can be shadowed for long periods of time. It is shown, in addition that the stretching of the phase space near escape and capture regions is a cause for the failure of the shadowing method. []

In a numerical study the authors of Kovács & Érdi (2009) compare the concept of permanent and transient chaos on the basis of the Sitnikov problem. They also present a link between the stickiness effect and the concept of chaotic scattering close to invariant tori. To this end the authors find a chaotic saddle of the Sitnikov problem and show its fractal structure. []

The authors of Bountis & Papadakis (2009) investigate the generalized Sitnikov problem where the third body moves perpendicular on a line through the center of mass of $N-1$ equally massed primaries that themselves move on a circle. The study includes a detailed investigation of the intervals of stability and instability, the investigation of periodic solutions that are found off the z-axis, and finally provide the solution of the problem for N tending to infinity. []

In Hagel (2009) an original method of advanced perturbation theory is tested for the first time on the extended Sitnikov problem, where the the primaries are moving on a circle, and the mass of the third body coincides with the mass of one of the primaries. The resulting system of coupled equations is solved for small oscillations of the third body using function iteration, and advanced perturbative methods. []

The stability of motion of the third body when the primaries are oblate and radiating is studied in Kalantonis et al. (2008). Using perturbation theory based on Floquet theory the authors compute the stability of motion, they provide critical orbits at which families of periodic orbits bifurcate. It is shown that the family of straight line oscillations only exists for identical primary bodies. []

In Soulis et al. (2008), the authors study the existence and stability of straight line periodic orbits in a generalized Sitnikov problem, where three equally massed primaries move on circular orbit. They provide a detailed study of the stability interval along the z-axis, and also investigate the bifurcation of 3D families of symmetric periodic orbits by means of the SALI method and suitably chosen Poincaré maps. []

The stability intervals of motion of the third body along the z-axis are discussed, and extended into the x,y space in Soulis et al. (2007). The study is based on the classical Sitnikov problem with $m=0$, but also for small non-zero mass $m$. The authors find that the stability regions that are extended into 3D space grow as $|z|$ grows. Furthermore, the authors find, that for increasing $m$, the stability intervals close to $z=0$ disappear first, namely before the stability intervals situated at larger values of $z$. []

In Perdios (2007) the author investigates the stability and extension of the family of straight line periodic orbits into 3D, in the circular Sitnikov problem. Several new critical orbits are found at which families of three dimensional periodic orbits of the same or double period bifurcate. The study also investigates the influence of nearly equal primaries. []

A stroboscopic map is used in Kovács & Érdi (2007) to investigate the extended phase space of the standard Sitnikov problem. The authors calculate the escape times, and investigate properties of the phase space close to the origin and far away from it. []

A complete parameter study of the classical Sitnikov problem can be found in Dvorak (2007). The author provides detailed surfaces of sections in the parameter space $e$, and initial condition space $z(0)$. A shrinking of the main island with increasing $e$, and the dynamics of the 2-1 periodic orbit is demonstrated: the corresponding island disappears (reappears) by means of (inverse) pitchfork bifurcations. Furthermore, the role of sticky orbits, and escape channels on the dynamics are discussed. []

The linear stability of the family of periodic orbits, and critical orbits at which families of 3D periodic orbits bifurcate are presented in Perdios & Kalantonis (2006). The authors perform their study in the case of the photogravitational circular restricted three-body problem. []

An approximate solution of the elliptic Sitnikov problem is given also in Hagel & Lhotka (2005). The solution is obtained using Floquet theory, the Courant-Snyder transformation, and Poincaré-Lindstedt series expansions. The solution is valid for moderate eccentricities and initial amplitudes $z(0)$, and velocities $z'(0)$. The expansion order for the linearized solution is $O(e)^{17}$, the expansion order for the full nonlinear problem is $O(e,z)^7$. []

An approximate analytic solution that is valid for small bounded oscillations is derived in Faruque (2003). It is valid for small oscillations and moderate eccentricities of the primaries. The solution is obtained using Floquet theory, and the Courant-Snyder transformation. []

The Poincaré-Lindstedt perturbation method is used to obtain the solution of the circular problem in Faruque (2002). The analytic solutions are of 2nd order in the small parameters. []

In Jiménez-Lara Escalona-Buendía (2001) the authors calculate the symmetry lines of the Sitnikov problem and their dependency on the parameters using stroboscopic maps. They find families of periodic orbits and their bifurcations. Bifurcation diagrams and patterns of bifurcations are found too. []

Corbera & Llibre (2000) prove the existence of symmetric periodic orbits in the elliptic problem, and prove that not all the periodic orbits of the Sitnikov problem are symmetric orbits. The results rely on the presence of a Bernoulli shift as a subsystem of the Poincaré mapping of the Sitnikov problem. []

A global analysis of the generalized Sitnikov problem is done in Chesley (1999A). The author reduces the dynamical problem to a two-dimensional symplectic map that depends on the mass ratio of the primaries, and the angular momentum. A thorough investigation of the parameter space. and a representation in terms of symbolic dynamics can be found. The paper describes the regions of allowable and bounded motion, and discusses the dominant periodic orbit. []

The extended Sitnikov problem is investigated in Dvorak & Sun (1997), where the authors study the case of three equally-massed bodies that stay in the standard Sitnikov configuration, but such that the primaries are not moving on classical Keplerian orbits. The system is explored using qualitative arguments, and the aid of properly chosen surfaces of section. The authors find a stable 1:2 resonant orbit that bifurcates into unstable orbits with increasing energies. []

Regular solutions near the 3:2 resonant orbit are derived in Jalali & Pourtakdoust (1997) using a rotating coordinate system and the averaging method. Approximate solutions are found by the authors by means of Jacobian elliptic functions. The authors also provide a study of the breakdown of regular motion due to chaos for certain values of the eccentricity of the primary bodies. []

The work of Hevia & Ranada (1996) is interesting from an educational point of view: the authors use the Sitnikov problem to demonstrate that at certain values of the parameters chaos appears. They argue that it provides a very clear example for introducing students to dynamical systems that omit chaos. []

Analytical approximations for the Sitnikov problem are obtained and described in Wodnar (1995).[]

The solution of the circular Sitnikov problem, and its period, in terms of elliptic functions can be found in Belbruno et al. (1994). The authors study the linear stability of periodic orbits, and of the families of periodic orbits of the 3D circular restricted three body problem that bifurcate from them. The work includes a study of some bifurcated families close to the 1/2 mass ratio. []

A systematic numerical study of the dynamical problem can be found in Dvorak (1993). The author provides in this original article a parameter study, and shows the existence of invariant curves that exist up to a certain value of the initial conditions. []

In Martinez Alfaro & Chiralt (1993) the authors prove the existence of KAM invariant rotational curves for small $e$, and also show for larger eccentricities the existence of two complementary sequences of intervals of values of $e$ that accumulate to $1$, such that for one of the sequences KAM-tori around a fixed point persist. []

The results of a thorough reexamination of the original article of Sitnikov can be found in Wodnar (1993). The author corrects some of the original statements about the existence of oscillatory motions in the three body problem made by Sitnikov himself. []

A mapping, and a first integral of motion for the Sitnikov problem is derived in Hagel & Trenkler (1993) by making use of computer algebraic methods, that were originally designed for high energy particle accelerators. The results are based on the formulation of the problem in terms of an approximate polynomial differential equation, using Chebyshev approximation techniques. []

An analytic approach is used to obtain approximate analytic solutions of the Sitnikov problem in Hagel (1992). The author uses a transformation that reduces the linear part of the equation of motion to the type of a harmonic oscillator. In this setting it is possible to derive an approximate integral of motion that is then used to explicitly find the solution of the problem. []

A numerical study, that is based on Poincaré surface of sections and Lyapunov characteristic exponents is provided in Liu et al. (1991), Liu et al. (1991b). Here, the authors support a relation between the chaotic region and the eccentricity of the primary's orbit. []

A conformal mapping is derived and used in Liu & Sun (1990). On its basis, the authors are able to show the existence of a hyperbolic invariant set. They also measure the stochasticity of the mapping in terms of Kustaanheimo-Stiefel entropy. Good agreement is demonstrated between the analytical determinations and the results of numerical simulations. []

The continuation of the bifurcations into the case of unequal masses of the primaries, and some of the bifurcating families of the 3D periodic motions are computed and discussed in Perdios & Markelos (1987, 1988). []