Stability interchanges in a curved Sitnikov problem
Luis Franco-Pérez Marian Gidea Mark Levi Ernesto Pérez-Chavela


L.Franco-Pérez]Luis Franco-Pérez 
lfranco@correo.cua.uam.mx
M. Gidea]Marian  Gidea
marian.gidea@gmail.com
M. Levi]Mark Levi
levi@math.psu.edu
E. Pérez-Chavela]Ernesto  Pérez-Chavela
epc@xanum.uam.mx


[UAM]Departamento de Matemáticas Aplicadas y Sistemas, UAM-Cuajimalpa, Av. Vasco de Quiroga 4871,  México, D.F. 05348, México.
[YU]Department of Mathematical Sciences, Yeshiva University, 245 Lexington Ave, New York, NY 10016, USA.
[PSU]Mathematics Department, Penn State University, University Park, PA 16802, USA.
[ITAM]Departamento de Matemáticas, ITAM México, Río Hondo 1, Col. Progreso Tizapán, México D.F. 01080 .



We consider a curved Sitnikov problem, in which an infinitesimal particle moves on a circle under the gravitational influence of two equal masses in Keplerian motion within a plane perpendicular to that circle.  There are two equilibrium points, whose stability we are studying.  We show that one of the equilibrium points undergoes stability interchanges as the semi-major axis of the Keplerian ellipses approaches the diameter of that circle. To derive this result, we first formulate and prove a general theorem on stability interchanges, and then we apply it to our model.  The motivation for our model resides with the n-body problem in  spaces of constant curvature.





Control Allocation for Wide Area Coordinated Damping
    [
    A. Mahmoodzadeha.mahmoodzadeh@iau-boukan.ac.ir, B. MalekolkalamiB.Malakolkalami@uok.ac.ir 
==============================================================================================





§ INTRODUCTION





 §.§ A curved Sitnikov problem
 We consider the following curved Sitnikov problem:  Two bodies  of equal masses  (primaries)  move, under mutual gravity,  on Keplerian ellipses  about their center of mass. A third, massless particle is confined to a circle passing through the center of mass of the primaries, denoted by P_0, and   perpendicular to the  plane of motion of the primaries; the second intersection point of the circle with that plane  is denoted by P_1. We assume that the massless particle moves under the gravitational influence  of the primaries without affecting them. The dynamics of the massless particle has two equilibrium points, at P_0 and  P_1.  We  focus on the local dynamics near these two points, more precisely, on the dependence of the linear stability of these points on the parameters of the problem.

When the Keplerian ellipses are not too large or too small, P_0 is a local center and P_1 is a hyperbolic fixed point. When we increase the size of the Keplerian ellipses, as the distance  between P_1 and the closest ellipse  approaches zero, then P_1 undergoes stability interchanges.
That is, there exists a sequence of open, mutually disjoint intervals of values of  the semi-major axis of the Keplerian ellipses, such that, on each of these intervals the linearized stability of P_1 is strongly stable,  and each complementary interval contains values where the linearized stability  is not strongly stable, i.e., it is either hyperbolic or parabolic.
The length  of these intervals approaches zero when the semi-major axis of the Keplerian ellipses approaches the diameter of the circle on which the massless particle moves. This phenomenon is  the main focus in the paper.




It is stated in <cit.> and suggested by numerical evidence <cit.> that the
 linearized stability of the  point P_0 also undergoes stability interchanges when the size of the binary is kept fixed and the eccentricity  of the Keplerian ellipses approaches 1.
 

Stability interchanges of the type described above are ubiquitous in systems of varying parameters; they appear, for example,  in the classical Hill's equation and in the Mathieu equation <cit.>.  To prove the occurrence of this phenomenon in our curved Sitnikov problem, we first formulate a general result on stability interchanges for a    general class of simple mechanical systems. More precisely, we consider the motion of two bodies  — one massive and one massless — which are confined to a pair of curves and move under Newtonian gravity. We let the  distance between the two curves be controlled by some parameter λ.    We assume that the position of the infinitesimal particle that achieves the minimum distance between the curves is an equilibrium point. We show that, in the case when the  minimum distance between the two curves approaches zero, corresponding to λ→ 0, there existence a sequence of mutually disjoint open intervals  (λ_2n-1,λ_2n),  whose lengths approach  zero as λ→ 0, such that whenever λ∈ (λ_2n-1,λ_2n) the linearized stability of the equilibrium point is strongly stable,  and each complementary interval contains values of λ where the linearized stability  is not strongly stable.
From this result we derive the above mentioned stability interchange  result   for the curved Sitnikov problem.

The curved Sitnikov problem  considered in this paper is an extension of the classical Sitnikov problem described in Section <ref> (also, see e.g.,  <cit.>). When the radius of the circle approaches infinity, in the limit we obtain the classical Sitnikov problem — the infinitesimal mass moves along the line perpendicular to the plane of the primaries and passing through the center of mass. The equilibrium point P_1 becomes the point at infinity  and is of a degenerate  hyperbolic type. Thus,  stability interchanges of P_1 represent a new phenomenon that we encounter in the curved Sitnikov problem but not in the classical one.

Also in the last case, it is well known that for =0 the classical Sitnikov problem is integrable. In the case of the curved one, numerical evidence suggests that it is not  (see Figure <ref>).

The motivation for considering the curved Sitnikov problem  resides in  the n-body problem in spaces with constant curvature, and with models of planetary motions in binary star systems, as discussed in Section <ref>.



 §.§ Classical Sitnikov problem

We recall here the classical Sitnikov problem. Two bodies (primaries) of equal masses m_1=m_2=1 move in a plane on Keplerian ellipses  of eccentricity  about their center of mass, and a third, massless particle moves on a line perpendicular to the plane of the primaries and   passing through their center of mass. By choosing the plane of the primaries the xy-plane and the line on which the massless particle moves the z-axis, the equations of motion of the massless particle can be written, in appropriate units, as

    z̈=-2z/(z^2+r^2(t))^3/2,
 where r(t) is the distance from the primaries to their center of mass given by

    r(t)=1-cos u(t),

where u(t) is the eccentric anomaly in the Kepler problem. By normalizing the time we can assume that the period of the primaries is 2π, and

    r(t)=(1-cos t)+O(^2),

for small .

When   =0, i.e., the primaries move on  a circular orbit and the dynamics of the massless particle is described by a 1-degree of freedom Hamiltonian and so is integrable. Depending on the energy level, one has the following types of solutions:  an equilibrium solution, when the particle rests at the center of mass of the primaries; periodic solutions around the center of mass; escape orbits, either parabolic, that  reach infinity with zero velocity, or hyperbolic, that  reach infinity with positive velocity.

When ∈(0,1), the differential equation (<ref>) is non-autonomous and the system is non-integrable. Consider the case ≪ 1.
The system also has bounded and unbounded orbits, as well as  unbounded  oscillatory orbits and  capture orbits (oscillatory orbits are those for which   lim sup_t→±∞|z(t)|=+∞ and lim inf_t→±∞|z(t)|<+∞, and capture orbits  are those for which lim sup_t→ -∞|z(t)|=+∞ and
lim sup_t→ +∞|z(t)|<+∞).  In his famous paper about the final evolutions in the three body problem, Chazy introduced the term   oscillatory motions <cit.>, although he did not find examples of these, leaving the question of their existence open. Sitnikov's model yielded the first example of oscillatory motions <cit.>.
There are many relevant works on this problem, including <cit.>.

The curved Sitnikov problem introduced in Section <ref> is a modification of the classical problem when the  massless particle moves on a circle rather than a line.  Here  we regard  the circle  as  a very simple restricted model of a space with constant curvature. In Subsection <ref>
we introduce and summarize some aspects of  this problem.





 §.§ The n-body problem in spaces with

constant curvature

The   n-body problem on spaces with

constant curvature is a natural extension of the  n-body problem  in the Euclidean space; in either case  the gravitational law considered is Newtonian.  The extension was first proposed independently by the founders of hyperbolic geometry,  Nikolai Lobachevsky and János Bolyai. It was subsequently studied in the late 19th, early 20th century,  by Serret, Killing, Lipschitz, Liebmann, Schering, etc. Schrödinger developed a quantum mechanical analogue of the Kepler problem on the two-sphere in 1940. The interest  in the problem was revived  by     Kozlov, Harin,  Borisov,  Mamaev,  Kilin, Shchepetilov, Vozmischeva, and others, in the 1990's. A  more recent surge of interest  was stimulated by the works on relative equilibria in spaces with constant curvature (both positive and negative)  by Diacu, Pérez-Chavela, Santoprete,  and others, starting in the 2010's.  See <cit.> for a history of the  problem and a comprehensive list of references.

A distinctive aspect   of the n-body problem on curved spaces is that the lack of (Galilean) translational invariance results in the lack of  center-of-mass and linear-momentum integrals.
Hence, the study of the motion cannot be reduced to a  barycentric coordinate system.

As a consequence, the two-body problem on a sphere  can no longer be reduced to the
corresponding problem of motion in a central potential field, as is the case for   the Kepler problem in the Euclidean space. As it turns out, the two-body problem on the sphere is not  integrable <cit.>.

Studying the three-body problem on spaces with curvature is also challenging.
Perhaps  the simplest  model is the restricted  three-body problem  on a circle. This was studied in  <cit.>. First, they consider the motion of the two primaries on the circle, which  is integrable, collisions can be regularized, and all orbits can be classified into three different classes (elliptic, hyperbolic, parabolic). Then  they consider the motion of the massless particle  under the gravity of the primaries,   when one or both primaries are at a fixed position. They obtain once again a complete classification of all orbits of the massless particle.

In this paper we take the ideas from above one step further, by considering the curved Sitnikov problem,  with the  massless particle moving on a circle under the gravitational influence of two primaries that move on Keplerian ellipses in a plane perpendicular to that circle.   In the limit case, when the primaries are identified with one point, that is when the primaries coalesce into a single body, the Keplerian ellipses degenerate to a point, and  the limit problem coincides with the two-body problem on a circle described above.

While the motivation of this work is theoretical, there are possible connections with the dynamics of planets in binary star systems. About 20 planets outside of the Solar System have been confirmed to orbit about binary stars systems; since more than half of the main sequence stars have at least one stellar companion, it is expected that a substantial fraction of planets form in binary stars systems. The orbital dynamics of such planets can vary widely, with some planets orbiting one star and some others orbiting both stars. Some chaotic-like planetary  orbits have also been observed, e.g.  planet Kepler-413b orbiting Kepler-413 A and Kepler-413 B in the constellation Cygnus, which displays erratic precession. This planet's orbit is tilted away from the plane of binaries and deviating  from Kepler's laws. It is hypothesized that this tilt may be due to the gravitational influence of   a third star nearby <cit.>. Of a related interest is the relativistic version of the Sitnikov problem <cit.>.
Thus, mathematical models like the one considered in this paper could  be helpful to understand  possible types of planetary orbits in  binary stars systems.


To complete this introduction, the paper is organized as follows: In Section <ref> we go deeper in the description of the curved Sitnikov problem, studying the limit cases and its general properties. In  Section <ref> we present a general result on stability interchanges. In Section <ref> we show that the equilibrium points in the curved Sitnikov problem present stability interchanges. Finally, in order to have a self contained paper, we add an Appendix with general results (without proofs) from Floquet theory.



§ THE CURVED SITNIKOV PROBLEM



 §.§ Description of the model

We consider two bodies with equal masses (primaries) moving under mutual Newtonian gravity on identical elliptical orbits of eccentricity ε, about their center of mass. For small values of ε, the distance r(t) from either primary to the center of mass of the binary is given by

    r_ε(t;r)=    rρ(t;ε),      r>0,
    ρ(t;ε)=   (1-εcos (u(t)))= (1-εcos (t))+𝒪(ε^2),


where u(t) is the eccentric anomaly, which satisfies Kepler's equation u-sin u=(2π/τ)t, where τ/2π is the mean motion[The mean motion is the time-average angular velocity over an orbit.] of the primaries. The expansion in (<ref>) is convergent for < _c = 0.6627...; see <cit.>.

A massless particle is confined  on a circle of radius R passing through the center of mass of the binary and perpendicular to the plane of its motion. We assume that the only force acting on the infinitesimal particle is  the  component along the circle of the resultant of the gravitational forces exerted by the primaries. The motion of the primaries take place in the xy-coordinate plane and the circle with radius R is in the yz-coordinate plane. See Figure <ref>.

   

  We place the center of mass at the point (0,R,0) in the xyz-coordinate system. The position of the primaries are determined by the functions
  
    𝐱_1(t)   =   (r_ε(t;r)sin t ,R+r_ε(t;r)cos t ,0),
    𝐱_2(t)   =   (-r_ε(t;r)sin t ,R-r_ε(t;r)cos t ,0).


Note that t=0 corresponds to the passage of the primaries through the pericenter  at y=R± r(1-ε) and  t=π to the passage of the primaries through the apocenter at y=R± r(1+ε); both peri- and apo-centers lie on the plane of the circle of radius R in the y-axis.

  The position of the infinitesimal particle is 𝐱(t)=(0,y(t),z(t)) (taking into account the restriction of motion for the infinitesimal particle to the circle y^2+z^2=R^2). We will derive the equations of motion by computing the gravitational forces exerted by the primaries:
  
    𝐅_ε(y,t;R,r)   =   -𝐱-𝐱_1/||𝐱-𝐱_1||^3-
         𝐱-𝐱_2/||𝐱-𝐱_2||^3
       =   -(-r_ε(t;r)sin t ,y-R-r_ε(t;r)cos t ,z)/||𝐱-𝐱_1||^3
          -(r_ε(t;r)sin t ,y-R+r_ε(t;r)cos t ,z)/||𝐱-𝐱_2||^3 ,

  where the distance from the particle to each primary is
 
    ||𝐱-𝐱_1||=
         [r^2_ε(t;r)+2R^2+2Rr_(t;r)cos t -2y(R+r_ε(t;r)cos t )]^1/2,
    
         ||𝐱-𝐱_2||=[r^2_ε(t;r) +2R^2-2Rr_(t;r)cos t -2y(R-r_ε(t;r)cos t )]^1/2.



We note that when  r(1+ε)=2R the elliptical orbit of the primary with the  apo-center at y<R  crosses the circle of radius R, hence collisions between the primary and the infinitesimal mass are possible. Therefore we will restrict to r<2R/1+ε;  when ε=0, this means r<2R.


We write (<ref>) in polar coordinates, that is y=Rcos q, z=Rsin q, and we obtain
  
    𝐅_ε(q,t;R,r)   =   -(-r_ε(t;r)sin t ,Rcos q-R-r_ε(t;r)cos t ,Rsin q)/||𝐱-𝐱_1||^3
          -(r_ε(t;r)sin t ,Rcos q-R+r_ε(t;r)cos t  ,Rsin q)/||𝐱-𝐱_2||^3.


  The origin  q=0 corresponds to the point (0,R,0) in the xyz-coordinate system. Thus, the primaries move on elliptical orbits around this point.

  Next we will retain the component along the circle of the resulting force (<ref>). That is, we will ignore the constraint force  that confines the motion of the particle to the circle, as this force  acts perpendicularly to the tangential component of the gravitational attraction force. The unit tangent vector to the circle of radius R at (0,Rcos( q),Rsin( q)) pointing in the positive direction is given by 𝐮( q)=(0,-sin( q),cos( q)). The component of the force 𝐅_ε( q,t;R,r) along the circle is computed as
  
    𝐅_ε( q,t;R,r)·𝐮( q)=-(R+r_ε(t;r)cos(t))sin( q)/||𝐱-𝐱_1||^3-(R-r_ε(t;r)cos(t))sin( q)/||𝐱-𝐱_2||^3.


  The motion of the particle, as a Hamiltonian system of one-and-a-half degrees of freedom, corresponds to
  
    q̇   =   p ,
    ṗ   =   f_ε(q,t;R,r),

  where
  
    f_ε( q,t;R,r)   :=   𝐅_ε( q,t;R,r)·𝐮( q) ,
    
         ||𝐱-𝐱_1||   =   [ r^2_ε(t;r) +2R(1-cos  q)(R+r_ε(t;r)cos t) ]^1/2 ,
    
         ||𝐱-𝐱_2||   =   [ r^2_ε(t;r)
         +2R(1-cos  q)(R-r_ε(t;r)cos t)]^1/2 ,


Hence

    H_(q,p,t;R,r)=p^2/2+V_(q,t;R,r),

where the potential is given by  
    V_(q,t;R,r)   = -1/R(1/||𝐱-𝐱_1||+1/||𝐱-𝐱_2||)
    .






 §.§ Limit cases


The curved Sitnikov problem can be viewed as  a link between the classical Sitnikov problem and the Kepler problem on the circle, mentioned in  Section <ref>.



  §.§.§ The limit R→∞.


   We express (<ref>)  in terms of the arc length w=R q, obtaining

  
    f_ε(w,t;R,r)=-(R+r_ε(t;r)cos t )sin(w/R)/[r^2_ε(t;r) +2R(1-cos(w/R))(R+r_ε(t;r)cos t) ]^3/2
    
          -(R-r_ε(t;r)cos t )sin(w/R)/[r^2_ε(t;r) +2R(1-cos(w/R))(R-r_ε(t;r)cos t) ]^3/2

  which we can write   in a suitable form as
  
    f_ε(w,t;R,r)=-wsin(w/R)/(w/R)(1+r_ε(t;r)/Rcos t )/[r^2_ε(t;r) +2w^2(1-cos(w/R))/ ( w/R)^2(1+r_ε(t;r)/Rcos t)  ]^3/2
    
          -wsin(w/R)/(w/R)(1-r_ε(t;r)/Rcos t )/[r^2_ε(t;r) +2w^2(1-cos(w/R))/ ( w/R)^2(1-r_ε(t;r)/Rcos t)  ]^3/2.

  Letting R tend to infinity we obtain
   
    lim_R→∞f_ε(w,t;R,r)=-2w/(r^2_ε(t;r)+w^2)^3/2 ,

   which is the classical Sitnikov Problem.




  §.§.§ The limit r→0.


When we take the limit r→0 in (<ref>) we are fusing the primaries into a large mass at the center of mass and we obtain a  two-body problem on the circle.
The component force along the circle corresponds to

    lim_r→0f_ε( q,t;R,r)=-sin( q)/√(2)R^2(1-cos( q))^3/2 .


This problem was studied in  <cit.>   with a different  force given by

    -1/R q^2+1/R(2π- q)^2,

 the distance between the large mass and the particle is measured by the arc length (in that paper the authors assume that  R=1).
The potential of the force (<ref>) is

    V_1(q_1)=-1/R^2√(2)(1-cos (q_1))^1/2,

where q_1 denotes the angular coordinate, and the potential for (<ref>) is
      
    V_2(q_2)=-1/Rq_2-1/R(2π-q_2),

where q_2 denotes the angular coordinate.
Each problem  defines an autonomous system with Hamiltonian

    H_i(p_i,q_i)=1/2p_i^2+V_i(q_i) ,

taking p_1=dq_1/dt, p_2=dq_2/dt and i=1,2.
Let ϕ^i_t  be the   flow of the Hamiltonian  H_i, and let  A_i denote the phase space, i=1,   2


Using that all orbits are determined by the energy relations given by (<ref>), it is not difficult to define a homeomorphism g:A_1 → A_2 which maps orbits of system (<ref>) into orbits of system (<ref>). In the same way we can define a homeomorphism h:A_2 → A_1 which in fact is
g^-1. This shows the C^0 equivalence of the respective flows.
One can show that the two corresponding flows are C^0–equivalent.



We recall from <cit.> that the solutions of the two-body problem on the circle (apart from the  equilibrium  antipodal  to the fixed body)  are  classified in three families (elliptic, parabolic and hyperbolic solutions) according to their energy level. The elliptic solutions come out of a collision, stop instantaneously, and reverse their path back to the collision with the fixed body. The parabolic  solutions come out of a collision and approach the equilibrium as t →∞. Hyperbolic motions  comes out of a collision with the fixed body, traverse  the whole circle and return to a collision.

We remark that the two limit cases R→∞ and r→ 0 are not equivalent. Indeed,  in the case r→ 0 the resulting system is autonomous,
the point q=0 is a singularity for the system, and the  point q=π is a hyperbolic fixed point, while in the case R→∞  the resulting system is non-autonomous (for ≠ 0),    the point q=0 is a fixed point of elliptic type, and the point q=∞ is a degenerate hyperbolic periodic orbit.



 §.§ General properties




  §.§.§ Extended phase space, symmetries, and equilibrium points

It is clear that, besides the limit cases R→∞ and r→ 0, the dynamics of the system depends only on the ratio r/R, so we can fix R=1 and study the  dependence of the global dynamics on r where  0<r<2.
In this case using (<ref>)  and (<ref>) we get

    f_ε(q,t;r)=-(1+rρ(t;ε)cos(t))sin( q)/[r^2ρ^2(t;ε) +2(1-cos  q)(1+rρ(t;ε)cos t)
         ]^3/2
    
         -(1-rρ(t;ε)cos(t))sin( q)/[r^2ρ^2(t;ε) +2(1-cos  q)(1-rρ(t;ε)cos t)
         ]^3/2.



To study the non-autonomous system (<ref>) we will make the system autonomous by
introducing the time as an extra dependent variable

    𝒳_ε(q,p,s;r)={[         q̇          =          p;          ṗ          = f_ε(q,s;r);          ṡ          =          1 ]..



This vector field is defined on [0,2π]×ℝ×[0,2π], where we identify the boundary points of the closed intervals. The flow of 𝒳_ ε possesses symmetries defined by the functions
   
    𝕊_1(q,p,s)   =   (-q,-p,s) ,
    𝕊_2(q,p,s)   =   (q,p,s+2π) ,
    𝕊_3(q,p,s)   =   (q+2π,p,s) ,
    𝕊_4(q,p,s)   =   (q,-p,-s)

  in the sense that
         
         
  * 𝕊_1(𝒳_ε(q,p,s))=𝒳_ε(𝕊_1(q,p,s)),
         
  * 𝒳_ε(q,p,s)=𝒳_ε(𝕊_2(q,p,s)),
         
  * 𝒳_ε(q,p,s)=𝒳_ε(𝕊_3(q,p,s))
         
  * 𝕊_4(𝒳_ε(q,p,s))=-𝒳_ε(𝕊_4(q,p,s)),
      
  as can be verified by a direct computation.


The function 𝕊_1 describes  the  symmetry respect to the trajectory (0,0,s), 𝕊_2 and 𝕊_3  describe the bi-periodicity of  f_ε(q,s;r) and 𝕊_4 describes the reversibility of the system.

System (<ref>) has two equilibria  (0,0) and (π,0), which correspond to periodic orbits for 𝒳_ε.



While the classical Sitnikov equation     is autonomous for ε = 0, our equation (<ref>) is not, and thus we expect it to be non–integrable, as is borne out by numerical simulation.  Figure <ref>  shows a Poincaré section corresponding to s=0 (mod 2π), with ε=0 and r=1. This simulation suggests that the invariant KAM circles
coexist with chaotic regions.
   



In the sequel, we will  analyze the dynamics around the  equilibrium points (π,0) and (0,0).  One important phenomenon that we will observe is that both equilibrium points  undergo stability interchanges as parameters are varied. More precisely, when 
ε  sufficiently  small is kept fixed and r→2R/1+ε, the point (π,0) undergoes infinitely many changes in stability, and when r is kept fixed  and ε→ 1, the point (0,0) undergoes infinitely many changes in stability.

In the next section we will first prove a general result.



§ A GENERAL RESULT ON STABILITY INTERCHANGES


In this section we switch to a  more general mechanical model which exhibits stability interchanges. We consider  the motion under mutual gravity of an infinitesimal  particle  and a heavy mass each constrained to its own curve and moving under gravitational attraction, and  study the linear stability of the equilibrium point corresponding to the closest position between the particles along the curves they are moving on. In Section <ref>, we will apply this general result to  the equilibrium point P_1 of the curved Sitnikov problem described in Section <ref>.  The fact that in the curved Sitnikov problem there are two heavy masses, rather than a single one as considered in this section, does not change the validity of the stability interchanges result, since, as we shall see,  what it ultimately matters is the time-periodic gravitational potential acting on the infinitesimal  particle.  



To describe the setting of this section, consider   a particle   constrained to a curve x =  x (s, λ ) in ℝ  ^3, where s is the arc length along the curve and λ is a parameter with values in some interval [0, λ_0], λ_0> 0.
Another (much larger) gravitational mass undergoes a prescribed  periodic  motion according to  y =  y (t, λ) =   y (t+1, λ ); see Figure  <ref>. We assume that the mass of the particle at x (s) is negligible compared to the mass at  y(t), treating the particle at x (s) as  massless.



 

To write the equation of motion for the unknown coordinate s of the massless particle,  let:

    z (s,t, λ ) =  x (s, λ) -  y (t, λ)
    	.

Assume that s=0,   t=0 minimize the distance between the two curves:

    |  z (0,0, λ) | = min_s,t | z (s,t, λ)|def=δ( λ ),

for all λ∈ [0, λ_0],  and that this minimum point is non-degenerate with respect to t, in the sense that 
    ∂^2 /∂ t^2|z(0,t,λ)|_| t=0≠ 0 .


Moreover, we make the following orthogonality assumption:

    ẋ(0,λ)·ẏ(0,λ)=0,  for all λ∈[0,λ_0].

	


In the sequel we will study the case  when the minimum distance min_s,t | z (s,t, λ)| → 0, that is, the shortest distance from the orbit y(t,λ) of the massive body  to the curve x(s,λ)
drawn by the infinitesimal mass approaches 0 as λ→λ_0. See Figure <ref>.
In the curved Sitnikov problem, this corresponds to the case when r→2R/1+ε (r→ 2R when 
 ε  = 0).



To write the equation of motion for s, we note that
 the Newtonian gravitational potential of the particle at x (s) is a function of s and t  given by

    U( s, t, λ  )  =- |  z (s,t, λ )  | ^-1,

and the evolution of s is    governed by the Euler–Lagrange equation
d/dt L_ṡ-L_s=0 with the Lagrangian

    L = 1/2ṡ ^2 - U( s, t, λ  ),

leading to[to explain this form of the Lagrangian, we note that  the equations for our massless particle are   obtained by taking the limit of the particle of small mass m; for such a particle, in the ambient potential U,  the Lagrangian is
1/2 m ṡ ^2    -m U( s, t, λ  ) – the factor m in front of U is due to the fact that U is the potential energy of the unit mass. Dividing the Euler–Lagrange equation  by m gives   (<ref>).]

    s̈ + U ^' ( s , t,λ)  = 0,   where   ^' = ∂/∂ s .


Note that s=0 is an equilibrium for any λ, since U ^' (0,t, λ ) = 0 for all t and for all λ, according to   (<ref>).
Linearizing   (<ref>)
around the equilibrium s=0 we obtain

    S̈+ a(t, λ  ) S = 0,    a(t+1, λ ) = a (t, λ ),

where a(t, λ  ) =  U ^'' (0,t,λ).



We have the following general result:


 
Assume that  (<ref>), (<ref>),
 (<ref>)   hold, that x and y are both bounded in the C^2–norm uniformly in λ, and

    min_s,tx(s,λ)-y(t,λ)→ 0  as λ→ 0.

Then there exists an infinite sequence

    λ_1 > λ_2 ≥λ _3> λ  _4≥⋯ >
    	λ_ 2n-1>λ_ 2n≥⋯→ 0

such that  the equilibrium solution s=0 of   (<ref>)   is linearly  strongly stable[The equilibrium solution is linearly  strongly stable if the Floquet multipliers of   (<ref>)  lie on the unit circle and are not real, or equivalently, if the linearized system lies in the interior of the set of stable systems.] for all λ∈ (λ_2n, λ  _2n-1). Furthermore, each complementary λ–interval  contains points where the linearized equilibrium is not strongly stable, i.e. is either hyperbolic or  parabolic.



The proof of  Theorem <ref> relies on Lemmas <ref>, <ref> and <ref> stated below.

   Consider the linear system
 
    ẍ + a(t, λ ) x=0,

where a(t,λ) is a continuous function  of t∈ [0,1], λ∈ (0, λ̅], where λ̅> 0.  Let z(t, λ) = x+ i ẋ, where x=x(t, λ ) is    a nontrivial solution of   (<ref>). Assume that there exists an
  interval [t_0( λ ), t_1( λ ) ]⊂ [0,1], possibly depending on λ, such that

    lim_λ↓ 0 z(t, λ)|_t_0^t_1→ -∞.

Here z(t, λ ) is defined as a continuous function of t. Although this choice of  is unique only modulo 2 π,
its increment as stated in equation (<ref>) is uniquely defined.

Then there exists a sequence  {λ_k}_k=0^ ∞ satisfying (<ref>)  such that the Floquet matrix of   (<ref>)  is strongly stable
for all  λ∈ (λ_2n, λ  _2n-1), for any n > 0. Furthermore, every complementary λ–interval contains
values of λ for which the Floquet matrix is not strongly stable.


Proof of this lemma can be found in <cit.>.

  
Consider the linear system

    ẍ + a(t) x=0,
    and assume that a(t) > 0 on some interval
[t_0, t_1].  For any (nontrivial) solution x(t) the corresponding  phase vector
z(t ) = x + i ẋ rotates by

    θ [z]def= z(t) |_t_0^t_1≤ -min_[t_0,t_1] √( a(t) ) (t_1-t_0) + π.




Writing the differential equation ẍ + a (t) x = 0 as a system ẋ = y,  ẏ  = - a(t)x, we obtain  (using complex notation): 

    θ̇=  d/dt (x+iy)= - d/dtIm ( ln z) = - Imż/z =
    	-Im(y-iax)(x-iy)/x ^2 + y ^2  =  -(a cos^2  θ+ sin ^2 θ ).

We conclude that for any solution of   (<ref>), the angle   θ= (x+iy) satisfies

    θ̇≤  -(a_mcos ^2 θ + sin^2 θ ),   for    t∈ [t_0,t_1] ,   where    a_m=min_[t_0,t_1] a(t).

To invoke comparison estimates, consider θ̅(t) which satisfies

    d/dtθ̅=  -(a_m cos^2 θ̅+ sin^2 θ̅),   θ̅(t_0) = θ (t_0).

By the comparison estimate, we conclude:

    θ|_t_0^t_1≤θ̅|_t_0^t_1,

 and the proof of the lemma will be complete once we show that θ̅ satisfies the estimate   (<ref>).
To that end we  consider a solution of

    ẍ̅̈ + a_m x̅ = 0

with the initial condition satisfying

    (x̅(t_0)+iẋ̅̇(t_0)) = θ̅(t_0).
    This  solution is of the form

    x̅(t)= Acos  (√( a_m) t- φ ),    A = const. ,

where φ is chosen so as to satisfy   (<ref>).

Since ( x̅+ iẋ̅̇) satisfies
the same differential equation as θ̅, and since the  initial conditions match, we conclude that
θ̅(t)  = ( x̅+ iẋ̅̇), so that
    θ̅(t)  =
    	( cos  (√( a_m )  t- φ) - i√(  a_m )sin (√(  a_m ) t-  φ ) )=
    	- (√(a_m) t- φ )+ π /2,

where X denotes a quantity whose absolute value does not exceed X.  In other words, θ̅(t) is given by a linear function with coefficient -  √( a_m), up to an error
< π /2.
The last inequality is due to the fact that the complex numbers cos(√(a_m)t-ϕ)-isin(√(a_m)t-ϕ) and
cos(√(a_m)t-ϕ)-i√(a_m)tsin(√(a_m)t-ϕ) lie in the same quadrant, so the difference between their arguments is no more than π/2.


Therefore, over the interval [t_0,t_1] the function θ̅ changes by  the amount √( a_m)(t_1-t_0) with the error of at most π/2 + π/2 = π:

    θ̅|_t_0^t_1≤ -√(  a_m ) (t_1-t_0) + π;

restating this more formally,   (<ref>)  implies

    θ̅(t_1)<  -  (√( a_m) t_1- φ)+ π/2 ,    θ̅(t_0)>  -  (√( a_m)t_0- φ)- π/2.

Subtracting the second inequality from the first gives   (<ref>).

Substituting  (<ref>) into  (<ref>) yields (<ref>)  and completes the proof of  Lemma <ref>. 


 
Consider the potential U defined by equation (<ref>).    Assume that the functions
x= x(s,λ), y= y(t,λ) satisfy


  (i)   the minimum min_s,t | z(·,·,λ) |=δ(λ) is non-degenerate and   is achieved at s=t=0, where z(s,t,λ)= x(s,λ)- y(t,λ),


  (ii)  || z(·,·,λ) ||_C^2≤ M  uniformly in 
0<λ < λ_0, and


  (iii) δ ( λ ) → 0  as λ→ 0.


Then there exists time τ= τ (λ) (approaching zero as λ→ 0) such that


    lim_λ→ 0(τ(λ) ·min_|t|≤τ (λ)√(U ^'' (0,t,λ)))	→∞.




Differentiating   (<ref>)  with respect to s twice, we get

    U^'' (0,t, λ )   = [( z^'· z^' +  z· z  ^'')
    	( z· z )  -3( z· z ^' )^2 /( z· z)^5/2]_| s=0.

We now estimate all the dot products in the above expression to obtain a lower bound.

First,

    z ^'· z ^' = 1,

since  s is the arc length, and from here z ^'· z ^'' = 0. Now, to  estimate
z· z and z· z^' we observe that  z· z^'= 1/2 ( z· z)^'  and we note that the first expression, as a  function of t with s=0 fixed has a   minimum at t=0 that we call δ ^2, and that the second function vanishes at t=0.  Applying Taylor's formula with respect to t we then have
 
    z· z  =δ^2 + k  t ^2,  z· z^' =  1/2
    	( z· z)^'  =  kt ,

where the constant k is determined by  the C^2–norm M of z. For the remaining dot product we have (still keeping s=0 and t arbitrary):

    | z· z ^''| ≤ |  z | | z ^''  |      (<ref>) ≤  M
    	√(δ ^2 + k t ^2 ).




Using the above estimates in  (<ref>),  we obtain

    U ^''(0,t, λ ) ≥(1- M √(δ^2 + k     t ^2 ))δ^2- 3k ^2 t^2 /(δ^2 + k  t ^2)^5/2 .


Now we restrict  t to have  δ ^2 + kt ^2≤ 2 δ^2; this  guarantees that the denominator in  (<ref>) does not exceed (2 δ )^5/2; to bound the numerator, we further restrict t so that the dominant part
δ ^2 - 3 k ^2 t^2 ≥1/2δ ^2, thus bounding
the numerator from below by

    (  δ ^2 - 3 k ^2 t ^2 )  -  M  √(2 δ ^2)δ^2 ≥1/2δ ^2 - M  √(2)δ ^3 > 1/4δ ^2

if δ is sufficiently small.
Summarizing, we restricted t to

    | t | ≤ c δdef=τ ( λ ),   where  c = min(k^-1/2 ,
    	(k  √( 6) ) ^-1  ),

and showed that for all such t and for δ small enough

    U^''(0,t, λ ) ≥1/4δ ^2  /(2δ^2 )^5/2  = c_1/δ ^3,

where c_1=2 ^-9/2.
With τ defined in   (<ref>) we obtain
  lim_τ→ 0( τ·min_|t| ≤τ√(a(t,λ))) = ∞, thus completing the proof of the lemma.
 


Proof of Theorem <ref>
Consider  the linearized equation

    ẍ + U^'' (0,t, λ)x = 0,


and consider the phase point   z=x+i ẋ of a nontrivial solution.
 Lemma <ref> gives us the rotation estimate (<ref>) for any time interval [t_0, t_1]; let this interval be
 [- τ ( λ ) ,τ ( λ )] where τ ( λ )  is taken from the statement of Lemma  <ref>. We then have  from   (<ref>):
  (<ref>)
 	
    θ [z]def= z(t) |_t_0^t_1≤
    	- 2τ(λ )·min_|t|≤τ(λ )√(U^''(0,t, λ)) + π.

According to the conclusion  (<ref>)  of   Lemma <ref>,
θ→ - ∞ as λ→ 0. This satisfies the condition   (<ref>)
Lemma <ref>, which now applies and its conclusion comples the proof of  Theorem <ref>.  




§ STABILITY OF THE EQUILIBRIUM POINTS IN THE CURVED SITNIKOV PROBLEM

The system (<ref>) has two equilibria  (0,0) and (π,0), which correspond to periodic orbits for 𝒳_ε.
The associated linear system around the fixed point (q_*,p_*) can be
written as
   
    𝐯̇=A(t)𝐯 ,    𝐯=([ x; y;  ]) ,

   with
   
    A(t)=.(
                   [        0        1; ∂f_ε/∂ q        0;          ])|_q=q_*.


Let X(t) be a fundamental matrix solution of system (<ref>) given by
   
    X(t)=(
               [ x_1(t) x_2(t); y_1(t) y_2(t);        ]),

     with the initial condition X(0)=I, the identity matrix;  x_1 is an even function and x_2 and odd one since ∂f_ε/∂ q is an even function with respect to t. The monodromy matrix  is given by X(2π) and we denote λ_1,λ_2 its eigenvalues, the Floquet multipliers associated to (<ref>).
   These are given by
   
    λ_1,λ_2=x_1(2π)+y_2(2π)±√((x_1(2π)+y_2(2π))^2-4)/2 ,

   the trace  Tr(X(2π))=x_1(2π)+y_2(2π) determines the linearized dynamics around the fixed point. Moreover, since the function (∂f_ε/∂ q)_| q=q^*  is an even function, we know x_1(2π)=y_2(2π) (see <cit.>) and then
   
    λ_1,λ_2=y_2(2π)±√((y_2(2π))^2-1) .

   
   To emphasize the dependence on the parameters , r we write
   
    y_2(2π;ε,r)=y_2(2π).

 Thus the linear stability of the equilibrium is
   
      
  * Elliptic type: | y_2(2π;ε,r) |<1.
      
  * Parabolic type: |y_2(2π;ε,r)|=1.
      
  * Hyperbolic type: |y_2(2π;ε,r)|>1.
   
   Since the Wronskian is equal to 1 for all t, then
   
    W(ε, r)=(y_2(2π;ε,r))^2-x_2(2π;ε,r)y_1(2π;ε,r)=1,

   and we have:
   
    (y_2(2π;ε,r))^2=1+x_2(2π;ε,r)y_1(2π;ε,r) .

   From this last expression it follows  that, in  the parabolic case, the periodic orbit corresponding to the equilibrium point (π,0) is associated to   x_1(t) or to x_2(t) (or to both).





 §.§ Stability of the equilibrium point (π,0).  Case =0.
 In this subsection we consider the system (<ref>)-(<ref>) for the case ε=0, namely when the primaries are following circular trajectories. We also fix R=1. The extended vector field under study is
   
    𝒳_0(q,p,t;γ)={[         q̇          =          p;          ṗ          = f_0(q,t;r);          ṫ          =          1 ].

  where
  
    f_0(q,t;r):= -(1+rcos(t))sin(q)/||𝐱-𝐱_1||^3-(1-rcos(t))sin(q)/ ||𝐱-𝐱_2||^3 ,

     
    ||𝐱-𝐱_1||   =   [r^2+2+2 rcos(t)-2cos (q)-2rcos(q)cos(t))]^1/2 ,
    
         ||𝐱-𝐱_2||   =   [r^2+2-2 rcos(t)-2cos (q)+2rcos(q)cos(t))]^1/2 ,

  with 0<r<2. We observe that in this case, function (<ref>) is π–periodic (remember that f_ε(q,t;r) is 2π–periodic if ε >0).



   The linear system associated to (<ref>) around the fixed point (π,0) is defined by the function
   
    ∂f_0/∂ q(π,t;r)=(1+rcos(t))/[r^2+4+4rcos(t)]^3/2
         +(1-rcos(t))/[r^2+4-4rcos(t)]^3/2 ,

   which is C^1 respect to t and r.

 The function ∂f_0/∂ q(π,t;r)  is monotone decreasing with respect to   r, that is,  for all t∈[0,π)
   
    ∂f_0/∂ q(π,t;r_1)>∂f_0/∂ q(π,t;r_2)

   if r_1<r_2 (see Figure <ref>).


 Let be F(t,r)=∂f_0/∂ q(π,t;r), then
by straightforward computation we get

    ∂ F/∂ r(t,r)    =   (cos(t))[r^2+4rcos(t)+4] - 3(1+rcos(t))[r+2cos(t)]/[r^2+4+4rcos(t)]^5/2
       +   (- cos(t))[r^2-4rcos(t)+4] - 3(1-rcos(t))[r-2cos(t)]/[r^2+4-4rcos(t)]^5/2
       ≤   -4rcos(t)^2 - 6r/min{[r^2+4+4rcos(t)]^5/2,[r^2+4-4rcos(t)]^5/2} < 0.






   


   The equilibrium point (π,0) is of hyperbolic type if r≤(√(17)-3)^1/2=1.059….

   
We first show that F(t,r)≥ 0 for all t if and only if r≤(√(17)-3)^1/2.
We compute 
    ∂ F/∂ t(t,r)=rsin t [-1/(r^2+4rcos t+4)^1/2 +1/(r^2-4rcos t+4)^1/2.
    .
    +6(1+rcos t)/(r^2+4rcos t+4)^5/2-6(1-rcos t)/(r^2-4rcos t+4)^5/2 ].


We have F(t,r)=F( π-t, r ) so it is enough to restrict t∈[0,π/2]. Note that ∂ F/∂ t(t,r)=0 for t=0,π/2.

For 0<r<2, we have 
    F  (π/2,r )=2/(r^2+4)^3/2>0,

and 
    F(0,r)=1+r/(2+r)^3+1-r/(2-r)^3=-2(r^4+6r^2-8)/(4-r^2)^3.


It follows immediately that F(0,r)<0 if r>(√(17)-3)^1/2.
If rcos t≤ 1 then (<ref>) implies F(t,r)≥ 0. Hence F(t,r)≥ 0 for all t provided r≤ 1.
Let 1<r≤ (√(17)-3)^1/2. If t∈[0,cos^-1(1/r)], which is equivalent to rcos t>1, then 
    1/(r^2-4rcos t+4)^1/2
    -1/(r^2+4rcos t+4)^1/2≥ 0,
    6(1+rcos t)/(r^2+4rcos t+4)^5/2-6(1-rcos t)/(r^2-4rcos t+4)^5/2≥ 0.

Therefore (<ref>) implies that  F(t,r) is increasing in t for t∈ [0,cos^-1(1/r)], and, since F(0,r)≥ 0, it follows that F(t,r)≥ 0 for all t.

Now let  x(t)=x_1(t)  + x_2(t) be a particular solution of system (<ref>); by hypothesis it satisfies x(0)=1, ẋ(0)=1. From (<ref>) with ε = 0 we obtain

    ẍ = ∂f_0/∂ q(π,t;r) x    or  ẍ - F(t,r)x = 0.


Then, since  F(t,r)≤ 0 for all t, we can apply directly Lyapunov's instability criterion (see for instance page 60 in <cit.>) to show that (π,0) is of hyperbolic type.  

 We can in fact estimate the first value r_1 of r at which the equilibrium point (π,0) becomes of parabolic type for the first time. 
 The idea is to extend slightly the result beyond r=(√(17)-3)^1/2, and then find the maximum of the r for which Proposition <ref> holds. In the way we have to do straightforward analytic but tedious computations that we decide to avoid in this paper. Finally we get that the value of r to have the first parabolic solution is  r_1 ≈ 1.2472⋯.





   


In Figure <ref> we show a couple of numerical simulations which illustrate the stability interchanges of the equilibrium point (π,0) for r ∈ (0,2) and ε = 0. The figure on the right hand side is a plot for r close to 2.



 §.§ Stability of the equilibrium point (π,0).  Case ≠ 0.


In this case, for every 0<r<(√(17)-3)^1/2, ∂ f_0/∂ q (π,t;r)>0,  hence ∂ f_/∂ q (π,t;r)>0 for all >0 sufficiently small (depending on r). Therefore  (π,0)  remains an equilibrium point of hyperbolic type in the case when the primaries move on Keplerian ellipses of sufficiently small eccentricity ε>0.

We remark   that  Theorem <ref>  does not depend of the shape of the curves where x(x) or y(t) are moving, in other words we can apply Theorem <ref> independently if the primaries are moving on a circle or on ellipses of eccentricity ε>0. Thus we obtain the following result on  stability interchanges:






















In the curved Sitnikov problem, let us fix any ε∈ [0,1),   let R=1, and consider  r (the semi-major axes of the Keplerian ellipses traced out by the primaries) as the parameter. As r approaches 2/1+ε, the distance between a Keplerian ellipse and the circle of the massless particle approaches zero. There exists a sequence   r_n↑2/1+ε satisfying

    r_0 ≤ r_1 < r _2 ≤  r _3 ⋯ <
    	r_2n≤ r_2n+1 < ⋯,

such that the equilibrium point (π, 0) of equation (<ref>)
is strongly stable for r ∈ (r_2n-1, r _2n)
 and not strongly stable for some r  in the complementary intervals (r_2n , r _2n+1). In other words, the equilibrium point (π, 0) loses and then regains its strong stability infinitely many times as r increases towards 2/1+ε.



We verify that Theorem <ref> applies. The curve x(s,λ) of Theorem <ref> is represented by the circle of radius R=1,  the curve y (t,λ) is represented by the orbit of the primary that gets closer to (π,0) (when ε =0 the primaries are co-orbital), and the parameter λ corresponds to r.  The planes of  the two curves are perpendicular, as in (<ref>), and the minimum distance between the curves, given by (<ref>), is non-degenerate as in (<ref>), and it corresponds to the infinitesimal mass being at y=-1 and the closest primary to this point being at y=1-r(1+ε). Thus, the minimum distance is δ(r)=2-r(1+ε),  and it approaches 0 when r→2/1+ε. To apply Theorem <ref> we only need to verify that x and y are bounded in the C^2-norm uniformly in r. This is obviously true for y(t,r) since the motion of the primary is not affected by the motion of the infinitesimal mass; it is also true for x(s,r) since s=arclength and the motion lies on the circle of radius R=1, hence x(s,r)=1 (the radius of the circle), x^'(s,r)=1 (the unit speed of a curve parametrized by arc-length), and x^''(s,r)=1 (the curvature of the circle). Hence the conclusion of Theorem <ref> follows immediately.
 








 §.§ Stability of the equilibrium point q=(0,0).


The  linear stability of  (0,0) is the same as of the barycenter in the classical Sitnikov problem, since

    ∂f_ε∂ q(0)=-2R/r^3ρ(t;ε)^3.

The stability of this point can be treated very similarly to that of the origin for Hill's equation, so in this analysis we  use results from that theory.

As in the study of the other equilibrium point we start with the case  ε = 0. Here the function f_0(q,t;r)
 defined in equation (<ref>) around q=(0,0) is given by
  
    ∂ f_0/∂ q(q)=-2/r^3q+9+r^2+9r^2(cos(t))^2/3r^5q^3+𝒪(q^5) .

  The local dynamics is determined by the linear part. The eigenvalues are on the unit circle, given by σ_1,2=± i√(2/r^3), and the Floquet multipliers, which come from the monodromy matrix
  
    X(π)=
         (
           [             cos(√(2/r^3)π)   (√(r^3/2))sin(√(2/r^3)π); -(√(2/r^3)π)sin(√(2/r^3)π)              cos(√(2/^3)π);                            ]),      X(0)=I ,

  are λ_1,2=e^± i√(2/r^3)π.

Hence q=(0,0) is of elliptic type if √(2/r^3)≠ k for any k∈ℤ, and it is of parabolic type if √(2/r^3)=k for some k∈ℤ. In fact, in the parabolic case, if √(2/r^3)=2m, m∈ℤ, then  there exists a π-periodic solution, and if √(2/r^3)=2m+1, m∈ℤ, then there exists a 2π-periodic solution.



   The equilibium point q=(0,0) of  the system defined by   (<ref>)  is stable  for all r∈(0,2).

   
      The linear part possess Floquet multipliers which place the system in either the elliptic or the parabolic case.
      In the last case, when √(2/r^3)=k for some k∈ℤ, there exists two independent eigenvectors associated to the eigenvalues λ_1,2. Then the conjugacy  of
      A and ± I implies that the monodromy matrix is A = ± I and q=(0,0) is stable (see <cit.> for more details).

When the origin is of elliptic type for the linear system,
we observe that the coefficient of the term of order 3 in equation (<ref>) is greater than zero for all r∈(0,2), then  we can apply directly Ortega's theorem (see Appendix and <cit.>), and therefore we obtain that the equilibrium point (0,0) is stable for the whole system, that is, including the nonlinear part. 


We note that, in the case when   ε=0, stability interchanges in a weak sense appear as r→ 0,  since (0,0) switches between  elliptic type when r≠ (2/k^2)^1/3 and parabolic type when r= (2/k^2)^1/3, k∈ℤ^+, as noted before.



In the case when ε≠ 0, as we mentioned earlier, the linear stability of (0,0) is the same as in the classical Sitnikov problem, so it only depends on the eccentricity parameter ε.
The papers <cit.> state that there are stability interchanges when the size r of the binary is kept fixed and the eccentricity ε of the Keplerian ellipses approaches 1.

We should point out that Theorem <ref> does not apply to this case since the function  r_ε(t,r) describing the motion of the primaries — corresponding to y(t,λ) in Theorem <ref> —   does not remain bounded in
the C^2 norm uniformly in ε.










§ FLOQUET THEORY
 In order to have a self contained paper we add this appendix with the main results on Floquet theory, must of them are very well known for people in the field.

 Consider the linear system
  
    𝐱̇=A(t)𝐱 ,    𝐱∈ℝ^2 ,

  where A(t) is a T-periodic matrix-valued function. Let X(t) be the fundamental matrix solution
  
    X(t)=(
                [ x_1(t) x_2(t); y_1(t) y_2(t);        ])

  with the initial condition X(0)=I, the identity matrix. Let λ_1 and λ_2 be the eigenvalues (Floquet multipliers) of the monodromy matrix matrix X(T)   and let μ_1,μ_2 (Floquet exponents) be such that λ_1=e^μ_1T, λ_2=e^μ_2T.

  [Floquet's theorem]
     Suppose X(t) is a fundamental matrix solution for (<ref>), then
        
    X(t+T)=X(t)X(T)

     for all t∈ℝ. Also there exists a constant matrix B such that e^TB=X(T) and a T-periodic  matrix P(t), so that, for all t,
        
    X(t)=P(t)e^Bt.

  

  
     Let be λ a Floquet multiplier for (<ref>) and λ=e^μ T, then there exists a nontrivial solution x(t)=e^μ tp(t), with p(t) a T-periodic function. Moreover, x(t+T)=λ x(t).
  
  Thus, the Floquet multipliers lead to the following characterization:
  
     
  * If |λ|<1⇔Re(μ)<0 then x(t)→ 0 as t→∞.
     
  * If |λ|=1⇔Re(μ)=0 then x(t) is a pseudo-periodic, bounded solution. In particular when λ=1 then x(t) is T-periodic and when λ=-1 then x(t) is 2T-periodic.
     
  * If |λ|>1 then Re(μ)>0 and therefore x(t)→∞ as t
         →∞, an unbounded solution.
  

  
     If  the Floquet multipliers satisfy λ_1≠λ_2,  then the equation (<ref>) has two linearly independent solutions
         
    x_1(t)=p_1(t)e^μ_1t ,     x_2(t)=p_2(t)e^μ_2t ,

         where p_1(t) and p_2(t) are T-periodic functions and μ_1 and μ_2 are the respective Floquet exponents.
  
  In this way, the stability of the solution to (<ref>) is
  
     
  * Asymptotically stable if |λ_i|<1 for i=1,2.
     
  * Lyapunov stable if |λ_i|≤1 for i=1,2, or if |λ_i|=1 and  the algebraic multiplicity equals the geometric multiplicity.
     
  * Unstable if |λ_i|>1 for at least one i, or if |λ_i|=1 and the algebraic multiplicity is greater than the geometric multiplicity.
  

   A particular case for the equation (<ref>) is the so-called Hill's equation, namely the periodic linear second order differential equation:
   
    z̈+f(t)z=0 ,

   where f(t) is a π-periodic function. Equation (<ref>), as a first order system, is
  
    𝐯̇=A(t)𝐯 ,    𝐯=([ x; y;   ]) ,    
          A(t)=(
                   [     0     1; -f(t)     0;       ]) .

  Let (<ref>) be a fundamental matrix solution of (<ref>). The monodromy matrix corresponds to X(π) and, as before, let λ_1 and λ_2 be the Floquet multipliers associated to (<ref>) and μ_1 and μ_2 be the corresponding Floquet exponents.

  In this paper we assume that f(t) is an even function; from Floquet theory we know that x_1 is an even and x_2 is an odd function. Also, the trace of A(t) vanishes and
  
    λ_1λ_2=e^∫_0^πtr(A(t))dt=1 .

  Assuming μ_1=a+ib is the Floquet exponent  corresponding to the Floquet multiplier λ_1, the general solution to (<ref>) is characterized as follows:
  
     
  * Elliptic type: λ_1∈ℂ∖ℝ, with |λ_1|=1 (and λ_2=λ̅_1). The general solution is pseudo-periodic and can be written as
         
    𝐯(t)=c_1Re(𝐩(t)e^ibt/π)+c_2Im(𝐩(t)e^ibt/π) .

         The origin is Lyapunov stable.
     
  * Parabolic type: λ_1=λ_2=± 1.
        
           
  * If λ_1=1 and there are two linearly independent eigenvectors of the monodromy matrix, the general solution is
               
    𝐯(t)=c_1𝐩_1(t)+c_2𝐩_2(t) ,

               and is π-periodic and Lyapunov stable.
           
  * If λ_1=1 and there is just one eigenvector associated to this eigenvalue, thus the general solution is
               
    𝐯(t)=(c_1+c_2t)𝐩_1(t)+c_2𝐩_2(t) ,

               and is unstable.
           
  * If λ_2=-1 and there are two linearly independent eigenvectors of the monodromy matrix, the general solution is
               
    𝐯(t)=c_1𝐩_1(t)e^it+c_2𝐩_2(t)e^it ,

               and is 2π-periodic and Lyapunov stable.
           
  * If λ_1=-1 and there is just one eigenvector associated to this eigenvalue, thus the general solution is
               
    𝐯(t)=(c_1+c_2t)𝐩_1(t)e^it+c_2𝐩_2(t)e^it ,

               and is unstable.
        
     
  * Hyperbolic type: λ_1∈ℝ, but |λ_1|≠ 1 (and λ_2=1/λ_1).
         
            
  * If λ_1>1, then the solution is
            
    𝐯(t)=c_1𝐩_1(t)e^μ_1 t+c_2𝐩_2(t)e^-μ_1 t .

            
  * If λ_1<-1, then the solution is
            
    𝐯(t)=c_1𝐩_1(t)e^μ_1 te^it+c_2𝐩_2(t)e^-μ_1 te^it .

         
         Thus, the origin is unstable.
  




A useful result from Rafael Ortega <cit.>  considers the nonlinear Hill equation

    y”+a(t)y+c(t)y^2n-1+d(t,y)=0,

with n≥2,  where the functions a,c:ℝ→ℝ are continuous,  T-periodic, and
∫_0^T|c(t)|≠ 0,
and the function d:ℝ×(-ϵ,ϵ)→ℝ, for ϵ>0, is continuous, has continuous derivatives of all orders respect to y, is T-periodic function respect to t, and
d(t,y)=𝒪(|y|^2n)  as
y→0 uniformly with respect to t∈ℝ.

The linear part around the solution y=0 of (<ref>) is

    y”+a(t)y=0 .


[Ortega's Theorem]
   Assume the following:
   
      
  * The equation (<ref>) is stable.
      
  * c≥0 or c≤0.
   
   Then y=0 is a stable solution of (<ref>).




 §.§ Acknowledgements
 We thank to the anonymous referees, their remarks and suggestions help us to improve this paper. Research of M.G. was partially supported by NSF grant  DMS-1515851.
M. L. gratefully acknowledges support by the NSF grant DMS-1412542.
The fourth author (EPC) has received partial support by the Asociación Mexicana de Cultura A.C.
Parts of this work have been done while the authors visited CIMAT, Guanajuato, LFP and EPC visit Yeshiva University and M.G. visited UAM-I in Mexico City. All authors are grateful for the hospitality of these institutions.




99
Alekseev1968a V. Alekseev, Quasirandom dynamical systems I, Math. USSR Sbornik 5, 73-128 (1968).

Alekseev1968b V. Alekseev, Quasirandom dynamical systems II, Math. USSR Sbornik, 6, 505-560 (1968).

Alekseev1969 V. Alekseev, Quasirandom dynamical systems III, Math. USSR Sbornik, 7, 1-43 (1969).

    Ces  L. Cesari, Asymptotic behavior and instability problems in ordinary differential equations, Springer-Verlag, New York, 1971.

   Chazy J. Chazy, Sur l'allure finale du mouvement dans le problème des trois corps I, II, III. Ann. Sci. Ecole Norm. Sup. 3^e Sér. 39 (1922); J. Math. Pures Appl. 8 (1929); Bull. Astron. 8 (1932).

   Dankowicz H. Dankowicz, P. Holmes, The existence of transversal homoclinic points in the Sitnikov problem. J. Differential Equations 116, 468-483 (1995).

   Diacu F. Diacu, Relative Equilibria in the 3-Dimensional Curved-N-Body Problem, Memoirs of the American Mathematical Society, Vol. 228, 2014.

    Fr L. Franco-Pérez, E. Pérez Chavela, Global symplectic regularization for some restricted 3–body problems on S^1. Nonlinear Analysis, 71, 5131-5143 (2009).

   GarciaPerezChavela  A. Garcia, E. Pérez-Chavela, Heteroclinic phenomena in the Sitnikov problem, Hamiltonian systems and celestial mechanics (Patzcuaro, 1998), 174-185, World Sci. Monogr. Ser. Math., 6, World Sci. Publ., River Edge, NJ, 2000.

   GorodetskiKaloshin A. Gorodetski and  V. Kaloshin, Hausdorff dimension of oscillatory motions in the restricted planar circular three body problem and in Sitnikov problem, preprint.

   Hagel_Lothka J. Hagel, C. Lhotka, A High Order Perturbation Analysis
         of the Sitnikov Problem, Celestial Mechanics and Dynamical Astronomy 93, 201-228 (2005).

   Kalas V. O. Kalas, P. S. Krasil'nikov, On equilibrium stability in the Sitnikov problem,
   Cosmic Research 49,  534-537  (2011).

   Kostov2014 V.B. Kostov, P.R. McCullough, J.A. Carter, M. Deleuil, R.F. Díaz, D. C. Fabrycky, G. Hébrard, T.C. Hinse, T. Mazeh, J.A. Orosz, Z.I. Tsvetanov, W.F. Welsh, Kepler-413b: a slightly misaligned, Neptune-size transiting circumbinary planet,  Arxive: 1401.7275 (2014).

    Kovacs T. Kovács, Gy. Bene and T. Tél, Relativistic effects in the chaotic Sitnikov problem, Monthly Notices of the Royal Astronomical Society, 414  (3),   2275-2281  (2011).

    Levi M. Levi, Stability of the inverted pendulum - a topological explanation. SIAM Review, 30 (4), 639-644 (1988).

 Alfaro  J. Martínez Alfaro, C. Chiralt, Invariant rotational curves in Sitnikov's Problem,
Celestial Mechanics and Dynamical Astronomy, 55 (4), 351–367 (1993).

   McGehee R. McGehee, A stable manifold theorem for degenerate fixed points with applications to
    Celestial Mechanics, J. Differential Equations 14, 70-88 (1973).

   Mos J. K. Moser, Stable and Random Motion in Dynamical Systems. Annals Math. Studies 77, Princeton University Press, 1973.

    Ortega R. Ortega, The stability of the equilibrium of a nonlinear Hill's equation. SIAM J. Math. Anal. 25 (5), 1393-1401 (1994).

   Ortega2 R. Ortega, The stability of the equilibrium: a search for the right approximation. Ten Mathematical Essays on Aproximation in Analysis and Topology, (J. Ferrera, J. López-Gómez and F.R. Ruiz del Portal, Eds.), Elsevier, 215-234, (2005).

    Plummer  H.C. Plummer, An Introductory Treatise on Dynamical Astronomy, New York, Dover, 1960.

   Robinson C. Robinson, Homoclinic Orbits and Oscillation for the
Planar Three-Body Problem, Journal of Differential Equations 52, 356-377 (1984).

   Sitnikov K. Sitnikov, The Existence of Oscillatory Motions in the Three-Body Problem. Translation from Doklady Akademii Nauk SSSR,  133 (2), 647-650 (1961).

  Magnus M. Wilhelm, W. Stanley, Hill's Equation. First edition, Dover,  1979.

  Ziglin S.L. Ziglin, Non-integrability of the restricted two-body problem on a sphere. Doklady RAN. 379  (4),  477-478 (2001). Engl. transl.: Physics-Doklady.  46  (8),  570-571 (2001).


