Second-kind symmetric periodic orbits for planar perturbed Kepler problems and applications

Abstract

We investigate the existence of families of symmetric periodic solutions of second kind as continuation of the elliptical orbits of the two-dimensional Kepler problem for certain symmetric differentiable perturbations using Delaunay coordinates. More precisely, we characterize the sufficient conditions for its existence and its type of stability is studied. The estimate on the characteristic multipliers of the symmetric periodic solutions is the new contribution to the field of symmetric periodic solutions. In addition, we present some results about the relationship between our symmetric periodic solutions and those obtained by the averaging method for Hamiltonian systems. As applications of our main results, we get new families of periodic solutions for: the perturbed hydrogen atom with stark and quadratic Zeeman effect, for the anisotropic Seeligers two-body problem and to the planar generalized Størmer problem.

1. Introduction

The study of periodic solutions in celestial mechanics not only benefits the development of mathematics, but also provides the intermediate orbits for space missions [16]. In celestial mechanics, one of the most well-known integrable model is the Kepler problem. There exist many other problems that are formulated as a perturbation of the Kepler problem in Cartesian coordinates (see [5, 6, 10, 13] and references therein) or in rotating coordinates (see, e.g., [16, 21, 25, 27]). Poincaré [28] considered the investigation of periodic solutions of the restricted three-body problem where, in particular, he classified the periodic orbits of second kind that are generated by the elliptic orbits of the planar Kepler problem (the first kind are generated by the circular orbits of the planar Kepler problem). As in the restricted three-body problem, we call the solutions of the second kind in the perturbed Kepler problems those that are generated by the planar elliptic orbits of the Kepler problem. The objective of this paper is to show analytically the existence of several families of symmetric periodic solutions of second kind of differentiable systems which are symmetric perturbations of the planar Kepler problem in cartesian coordinates. In a future work, we intend to deal with the spatial case.

More precisely, we want to analyse when the second kind of periodic orbits of the planar Kepler problem can be or not extended to periodic orbits of symmetric perturbed Kepler problems whose Hamiltonian function has the form

(1.1)\begin{equation} H({\bf q},{\bf p},\epsilon)=H_{0}({\bf q},{\bf p})+ \epsilon^{\alpha} H_{1}({\bf q}, {\bf p})+ H_R({\bf q},{\bf p} ,\epsilon), \end{equation}

with ${\bf q}= (x, y) \in \mathbb {R}^2 \setminus \{(0, 0)\}$, ${\bf p}= (p_x, p_y) \in \mathbb {R}^2$, $\alpha \in \mathbb {N}$, where $H_{0}({\bf q},{\bf p})= \frac {\|{\bf p}\|^2}{2} - \frac {1}{\|{\bf q}\|}$ is the two-dimensional Kepler problem, and the perturbed functions $H_{1}({\bf q}, {\bf p})$ and $H_R({\bf q},{\bf p} ,\epsilon )$ are both differentiable and $H_R$ is of order ${\mathcal {O}}(\epsilon ^{\alpha +1})$. Furthermore, we show the approximation of the characteristic multipliers associated with the symmetric periodic solutions, so its type of stability is characterized.

The continuation method to find periodic solutions goes back to Poincaré who in [28] studied the existence of periodic solutions in the three-body problem using the method that now is called Poincaré's continuation method [32].

There are considerable works that have contributed to find the periodic solution for perturbations of an integrable problem (see [1, 7]). For a perturbed Kepler problem, some of these works apply the method of average of first kind (see, e.g., [3, 4, 14, 25, 27]). The technique of combining discrete symmetries of Hamiltonian and the Poincaré continuation method has been considered in other works for two-degree-of-freedom ($2$-DOF) and three-degree-of-freedom problems (see, e.g., [5, 10, 25, 31–35]).

With regard to the study of periodic orbits of the second kind in perturbed Kepler problems, we can cite [2, 25, 31] and references therein. There are standard results in celestial mechanics which state that every elliptical solution of period $T=2\pi p/q$ ($p$ and $q$ are relative prime positive integers) of the planar rotating Kepler problem with initial condition on the $x$-axis can be extended to the perturbed planar circular restricted $3$-body problem. These solutions are symmetric with respect to the $x$-axis and have period $\tau$ close to $T$ (see [2, 21] ). More generally, the elliptical $T$-symmetric (with respect to the $x$-and $y$-axis) periodic solutions, where $T$ is as in the previous sentence, and can be extended to any planar symmetric perturbation of the Kepler problem (see [5]). But, when we consider a symmetric perturbation of the Kepler problem in the fixed inertial coordinates, the continuation of the elliptical Keplerian solutions depends on the first approach term $H_1$ given in equation (1.1). So, we need to establish sufficient conditions to extend the second-kind periodic orbits and how to proceed in the verification of them in concrete problems.

To our knowledge, there is no analogous and systematic study on the existence of periodic solutions of the second kind for a planar Keplerian perturbation in an inertial frame using Delaunay variables. There are a significant number of authors who search for periodic solutions to perturbations of the planar Kepler problem (in inertial frame) using Delaunay coordinates, but the technique used is the average method. Important references in this regard can be found in [2, 3, 20, 25, 27]. Many other authors considered the study of symmetric periodic orbits of first or second kind for planar perturbations of the Kepler problem. For the first kind of symmetric periodic solutions, we recommend the reader the references [5, 10, 34, 35]. For the second kind of symmetric periodic solutions, we mention [2, 8, 21, 25, 33].

The study of symmetric periodic orbits of a perturbed Kepler problem (1.1) also can be performed using averaging theory of Hamiltonian systems combined with symplectic reduction (see reference [36] and, more specifically, theorem 2.6). In this work, the authors gave sufficient conditions for the continuation of symmetric Keplerian solutions to the complete perturbed problem from the point of view of averaging theory. The results obtained in [36] depend on the implicit function theorem of Arenstorf. As a consequence of this study, the results are valid only for some values of the perturbed parameter.

In this work, we combine the discrete symmetries of the Hamiltonian function and the Poincaré continuation method, using strongly the first approximation of the solutions of the full Hamiltonian system given by a variational system. The technique used in this research is similar to that used in [5] for the existence of first kind symmetric solutions of the perturbed problem (1.1). The main contribution of this paper is to provide sufficient conditions for the existence of families of second-kind symmetric periodic solutions for planar Keplerian perturbations in inertial frame and using Delaunay variables. In addition, we give information about the stability of these solutions.

Moreover, we analyse that certain elliptic Keplerian solutions, which can be continued by the analytic Poincaré's continuation method as in theorem 2.3), also can give periodic solution to the full problem by using the averaging method for Hamiltonian systems, i.e., Reeb's Theorem.

To carry out our results, we have organized the contents of the paper as follows. In § 2, we write the problem in Delaunay coordinates. Then, using the continuation method, we prove a theorem that gives us the sufficient conditions for the existence of two types of families of initial conditions such that it gives us second-kind symmetric periodic solutions for the Hamiltonian (2.1). The first type of initial conditions depends on one small parameter and gives us symmetric periodic solutions with the same period of the Keplerian orbit. The second type of initial conditions depends on two small parameters and the periodic solutions generated by them have period close (not necessarily fixed) to the elliptical Keplerian orbit. In order to complement our study, we finish § 2 analysing the symmetric periodic solutions given by theorem 2.3 and those periodic solutions obtained using the averaging method. We point out the relationship of the period solutions obtained from these two different techniques. Moreover, the nontrivial characteristic multipliers of the symmetric periodic solutions are characterized, so the type of stability can be deduced. In § 3, as an application of our results, we obtain new periodic solutions for the $2$-DOF generalized Størmer problem, hydrogen atom with Stark and quadratic Zeeman effect and anisotropic two-body problem under Seeliger's potential and we give some important information on these models. In fact, after checking the literature on the subject, we emphasize that our results are new for this kind of dynamics of these problems. Finally, in § 4, some conclusions and prospective work are mentioned.

2. Second-kind symmetric periodic solutions

In this section and following the classification of Poincaré [28], we investigate the existence of second-kind symmetric periodic solutions associated with the Hamiltonian system with Hamiltonian function as in (1.1).

During this work, we are going to assume that the Hamiltonian function $H$ in (1.1) is invariant under one or both of the following anti-symplectic reflections

\begin{align*} & S_1:(x,y,p_x,p_y)\longrightarrow ({-}x,y,p_x,-p_y),\\ & S_2:(x,y,p_x,p_y)\longrightarrow (x,-y,-p_x,p_y). \end{align*}

The fixed sets of the symmetries $S_1$ and $S_2$ are the Lagrangian subspaces $\mathcal {L}_1=\{(0,y,p_x,0) \in \mathbb {R}^4: y, p_x\in \mathbb {R}\}$ and $\mathcal {L}_2=\{(x,0,0,p_y) \in \mathbb {R}^4: x, p_y\in \mathbb {R}\}$, respectively. Note that if $\varphi (t,{\bf q},{\bf p})=(x(t),y(t),p_x(t),p_y(t))$ is a solution associated with the Hamiltonian (1.1), then $S_j\circ \varphi (t,{\bf q},{\bf p})$ is also a solution. In particular, if we consider an initial condition $({\bf q},{\bf p})\in \mathcal {L}_j$ such that $\varphi (T/2,{\bf q},{\bf p})\in \mathcal {L}_j$, then the solution $\varphi (t,{\bf q},{\bf p})$ is $T$-periodic and $S_j$-symmetric, that is, knowing two different points on the orbit, we are able to get periodic solutions.

To obtain our results, we use Delaunay variables (see more details, e.g., in [8, 21, 33]), because with these coordinates it is easy to characterize the solutions of the Kepler problem and they are well defined for elliptic solutions. The Delaunay planar variables are $(\ell,g)$ (angular variables modulus $2\pi$) and $(L,G)$ (radial variables), where

\[ L=\sqrt{a},\quad G^2=L(1-e^2), \]

$e$ ($0\leq e <1$) is the eccentricity of the Keplerian orbit and $a$ is the semi-major axis of the ellipse (see figure 1). The angular variable $\ell$ is the mean anomaly, $g$ is the argument of the perigee measured from the ascending node, $L$ is the action related with the semi-major axis $a$ and $G$ is the angular momentum. We point out that the domain of the Delaunay variables is the set $D=\{(\ell,g,L,G): \, L>0, l,g \in [0,2\pi ), 0<|G|< L\}$, and of course, it excludes the rectilinear motions and circular solutions and they are well defined in a neighbourhood of an elliptic orbit of the Kepler problem.

Figure 1. Schematic of angular Delaunay variables. The angle $g$ is the argument of perigee and the angle $\ell$ is the ratio of the sector swept to the total area of the ellipse normalized by $2\pi$, i.e., if $S$ is the area of the sector swept out from perigee and $A$ is the area of ellipse then $\ell =2\pi S/A$.

Now express (1.1) in mixed polar-nodal coordinates and after we change polar-nodal to Delaunay variables. See [9, 24] for a complete description of the process of changing Cartesian coordinates to Delaunay coordinates. In Delaunay variables, the Hamiltonian (1.1) takes the form

(2.1)\begin{equation} {\mathcal{H}}(\ell,g,L,G,\epsilon)={\mathcal{H}}_{0}(L)+ \epsilon^{\alpha} {\mathcal{H}}_{1}(\ell,g,L,G)+{\mathcal{H}}_R(\ell,g,L,G,\epsilon), \end{equation}

where ${\mathcal {H}}_R(\ell,g,L,G,\epsilon )={\mathcal {O}}(\epsilon ^{\alpha +1})$, and the Hamiltonian of the Kepler problem is as follows:

(2.2)\begin{equation} {\mathcal{H}}_0(L)={-}\frac{1}{2L^2}. \end{equation}

The perturbed Hamiltonian system associated with the Hamiltonian (2.1) is written as

(2.3)\begin{equation} \begin{aligned} \dot \ell & = {\mathcal{H}}_{ L}= \frac{1}{L^3}+ \epsilon^{\alpha} \frac{\partial\mathcal{H}_{1}}{\partial L}+{\mathcal{O}}(\epsilon^{\alpha+1}), \quad \dot{L}={-}{\mathcal{H}}_{\ell}={-}\epsilon^{\alpha} \frac{\partial\mathcal{H}_{1}}{\partial l}+{\mathcal{O}}(\epsilon^{\alpha+1}),\\ \dot{g} & = {\mathcal{H}}_{ G}=\epsilon^{\alpha} \frac{\partial\mathcal{H}_{1}}{\partial G}+{\mathcal{O}}(\epsilon^{\alpha+1}), \quad \dot{G}={-}{\mathcal{H}}_{ g}={-}\epsilon^{\alpha} \frac{\partial\mathcal{H}_{1}}{\partial g}+{\mathcal{O}}(\epsilon^{\alpha+1}). \end{aligned} \end{equation}

We will denote by $\varphi (t)=\varphi (t, {\bf Y}; \epsilon )= (\ell (t,{\bf Y};\epsilon ), g(t,{\bf Y};\epsilon ), L(t,{\bf Y};\epsilon ),G(t,{\bf Y};\epsilon ))$, the flow of the Hamiltonian system associated with (2.3) with initial condition ${\bf Y}= (\ell _0,g_0,L_0,G_0)$ and we propose the following approximation of the solution:

(2.4)\begin{equation} \begin{aligned} & \ell(t,{\bf Y};\epsilon)= \ell^{(0)}(t,{\bf Y}) + \epsilon^{\alpha} \ell^{(1)}(t,{\bf Y}) + {\mathcal{O}}(\epsilon^{\alpha+1}),\\ & g(t,{\bf Y};\epsilon)= g^{(0)}(t,{\bf Y}) + \epsilon^{\alpha} g^{(1)}(t,{\bf Y}) + {\mathcal{O}}(\epsilon^{\alpha+1}),\\ & L(t,{\bf Y};\epsilon)= L^{(0)}(t,{\bf Y}) + \epsilon^{\alpha}L^{(1)}(t,{\bf Y}) + {\mathcal{O}}(\epsilon^{\alpha+1}),\\ & G(t,{\bf Y};\epsilon)= G^{(0)}(t,{\bf Y}) + \epsilon^{\alpha}G^{(1)}(t,{\bf Y}) + {\mathcal{O}}(\epsilon^{\alpha+1}). \end{aligned} \end{equation}

Obviously, the solution of the unperturbed system represents the first approximation and it is given by

(2.5)\begin{equation} \varphi_{0}(t, {\bf Y})=(\ell^{(0)}(t,{\bf Y}), g^{(0)}(t,{\bf Y}),L^{(0)}(t,{\bf Y}),G^{(0)}(t,{\bf Y}))= \left( \frac{t}{L_0^3} + \ell_0, g_0, L_0, G_0\right). \end{equation}

While the second approximation is characterized by

(2.6)\begin{equation} \begin{aligned} & \ell^{(1)}(t,{\bf Y})=\displaystyle\int_{0}^{t}\frac{\partial {\mathcal{H}}_1}{\partial L}(\varphi_{0}(\tau,{\bf Y}))\,{\rm d}\tau,\quad L^{(1)}(t,{\bf Y})={-}\displaystyle\int_{0}^{t}\frac{\partial {\mathcal{H}}_1}{\partial \ell}(\varphi_{0}(\tau,{\bf Y}))\,{\rm d}\tau,\\ & g^{(1)}(t,{\bf Y})= \displaystyle\int_{0}^{t}\frac{\partial {\mathcal{H}}_1}{\partial G}(\varphi_{0}(\tau,{\bf Y}))\,{\rm d}\tau,\quad G^{(1)}(t,{\bf Y})={-}\displaystyle\int_{0}^{t}\frac{\partial {\mathcal{H}}_1}{\partial g}(\varphi_{0}(\tau,{\bf Y}))\,{\rm d}\tau. \end{aligned} \end{equation}

Lemma 2.1 Let $\varphi (t)$ be a solution of system (2.3), then the following characterization holds:

1. (1) $\varphi (t)$ hits the fixed set $\mathcal {L}_1$ $($associated with the symmetry $S_1)$ at time $t=T$ if $\varphi (T)$ is on the perigee or apogee (i.e., $\ell = 0\mod \pi$) and furthermore, the perigee is on the $y$-axis (i.e., $g=\pi /2 \mod \pi$).

2. (2) $\varphi (t)$ hits the fixed set $\mathcal {L}_2$ $($associated with the symmetry $S_2)$ at time $t=T$ if $\varphi (T)$ is on the perigee or apogee (i.e., $\ell = 0\mod \pi$) and furthermore, the perigee is on the $x$-axis (i.e., $g=0\mod \pi$).

Proof. The proof follows directly from the definition of the sets $\mathcal {L}_1$ and $\mathcal {L}_2$ and the geometrical interpretation of the Delaunay variables. See figures 1 and 2.

Figure 2. Representation of a $S_1$-symmetric elliptical orbit and the possible positions of the perigee (depending on $k$). Left: $k=0$ ($\ell =\pi$ and $g=\pi /2$). Right: $k=1$ ($\ell =0$ and $g=\pi /2$).

Next, we take an elliptic solution $\varphi _{0}(t,{\bf Y}_{0}^{j,k})$ of the Kepler problem (2.2) in Delaunay variables with the initial condition

(2.7)\begin{equation} {\bf Y}_{0}^{j,k}= (Y_1^{(0)}, Y_2^{(0)}, Y_3^{(0)}, Y_4^{(0)})\equiv (0, g_{0}^{j,k}, L_0, G_0)\in \mathcal{L}_{j}, \end{equation}

where $|G_0|< L_0$, $j=1,2$ and

\[ g_{0}^{1,k}= (2k+1)\pi/2, \quad g_{0}^{2,k}=k\pi,\quad k=0,1. \]

Note that the index $k$ determines the position of the perigee on the $x$-axis or $y$-axis (see figure 2).

It is clear that the solution $\varphi _{0}(t, {\bf Y}_{0}^{j,k})=\left (\frac {t}{L_{0}^3}, g_{0}^{j,k}, L_0, G_0\right )$ of the Kepler problem is elliptic, $T=2\pi L_0^3$-periodic and its initial condition ${\bf Y}_{0}^{j,k} \in L_j$ for each $j=1,2$. After that, we take a small and convenient perturbation (in the ‘directions’ of $L$ and $G$) of the initial condition ${\bf Y}_{0}^{j,k}$ of the previous unperturbed elliptical solution in the form

(2.8)\begin{equation} {\bf Y}^{j,k}= (0, g_{0}^{j,k}, L_0+\delta L,G_0+\delta G) \in L_{j}. \end{equation}

Of course, the solution of the Kepler problem with this initial condition, $\varphi _{0}(t,{\bf Y}^{j,k})$, is also elliptic and it is given by

(2.9)\begin{equation} \begin{aligned} l(t) & =(L_0+\delta L)^{{-}3}t, \quad L(t)=L_0+\delta L,\\ g(t) & = g_{0}^{j,k},\quad G(t)= G_0+\delta G,\\ \end{aligned} \end{equation}

with $|\delta G|\leq \delta L$.

Remark 2.2 Here we discuss the effect of the perturbation of the Delaunay variables $L$ and $G$ in the initial condition ${\bf Y}_{0}^{j,k}$, ${\bf Y}^{j,k}$ in terms of the position and the velocity in cartesian coordinates called $({\bf q}_0, {\bf p}_0)$ and $({\bf q}(0), {\bf p}(0))$, respectively. Initially, note that the elliptic Keplerian orbit $\varphi _{0}(t,{\bf Y}^{j,k})$ has semi-major axis $(L_0+\delta L)^2$. So, one increment in $\delta L$ gives us a perturbation on the initial position of the orbit $\varphi _{0}(t,{\bf Y}_0^{j,k})$. On the other hand, since ${\bf Y}^{j,k} \in \mathcal {L}_j$, it follows that the angular momentum of the elliptic orbit $\varphi _{0}(t,{\bf Y}^{j,k})$ satisfies $|G_0+\delta G|= \|{\bf q}(0) \times {\bf p}(0)\|= \|{\bf q}(0)\| \|{\bf p}(0)\|$ because ${\bf q}(0) \perp {\bf p}(0)$. Note that by the previous discussion, we have ${\bf q}(0)={\bf q}_0+\delta {\bf q}$. So, we get $|G_0+\delta G|=\|{\bf q}_0+\delta {\bf q}\|\|{\bf p}(0)\|$, and we point out that an increment in $\delta G$ gives us a perturbation on the initial velocity of the elliptical orbit $\varphi _{0}(t,{\bf Y}_0^{j,k})$.

Before stating our main result, we introduce some notation in order to simplify the computations. Let $X= (X_1, X_2)= (\delta g,\delta G)$ (independent variables) and $\varphi _0(t, (Y_1^{(0)}, Y_2^{(0)}+X_1, Y_3^{(0)}, Y_4^{(0)}+ X_2))$ be a solution of the Kepler problem with initial condition $(Y_1^{(0)}, Y_2^{(0)}+X_1, Y_3^{(0)}, Y_4^{(0)}+ X_2)$, and

(2.10)\begin{equation} \bar{H}(X)= \int_0^T {\mathcal{H}}_1(\varphi_{0}(t, (Y_1^{(0)}, Y_2^{(0)}+X_1, Y_3^{(0)}, Y_4^{(0)}+ X_2)))\,{\rm d}t, \end{equation}

and consider the matrix

(2.11)\begin{equation} A= \mathbb{J}\left( \frac{\partial^2 \bar{H}}{\partial X_i \partial X_j}\right)_{X=0}, \end{equation}

where $\mathbb {J}= \left (\begin {array}{cc} 0 & I \cr -I & 0 \end {array}\right )$ denotes the standard skew-symmetric matrix. Now, we are ready to state our main result which provides sufficient conditions for the existence of $S_j$-symmetric periodic solutions ($j=1$ or $j=2$) of Hamiltonian (1.1) as continuation of the elliptic solutions of the Kepler problem.

Theorem 2.3 Fix the energy level ${\mathcal {H}}_{0}=- \frac {1}{2 L_0^2}$ and the period $T=2\pi L_0^3$ of the elliptic Kepler solution. Suppose that the Hamiltonian function $H$ in (1.1) is $S_j$-symmetric for $j\in \{1, 2\}$. Let $\varphi _{0}(t, {\bf Y}^{j,k})$ be an elliptical solution of the Kepler problem as in (2.9). Assume that the following two conditions are satisfied:

(2.12)\begin{equation} \begin{aligned} & (a)\quad \int_{0}^{T/2}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}\,{\rm d}\tau\\ & \quad = \int_{0}^{T/2}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf Y}^{j,k}=(0, g_{0}^{j,k}, L_0, G_0)}\,{\rm d}\tau=0, \\ & (b) \quad\int_{0}^{T/2}\frac{\partial^{2} {\mathcal{H}}_{1}}{\partial G \partial \delta G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{{\bf Y}^{j,k}}={\bf Y}_{0}^{j,k}}\,{\rm d}\tau\\ & \quad = {\int_{0}^{T/2}}\frac{\partial^{2} {\mathcal{H}}_{1}}{\partial G \partial \delta G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf Y}^{j,k}=(0, g_{0}^{j,k}, L_0, G_0)}\,{\rm d}\tau\neq 0. \end{aligned} \end{equation}

(i)   Then for $\epsilon$ sufficiently small there are $1$-parameter families of initial conditions

\[ {\bf Y}_{\epsilon}^{j,k}={\bf Y}_{0}^{j,k}+{\bf Y}^{j,k}(\epsilon) , \quad {\bf Y}^{j,k}(\epsilon)=(0,0,\delta L^{j,k}(\epsilon), \delta G^{j,k}(\epsilon)), \]

parametrized by $\epsilon$ such that $\varphi (t,{\bf Y}_{\epsilon }^{j,k};\epsilon )=\varphi _{0}(t, {\bf Y}_{0}^{j,k})+{\mathcal {O}}(\epsilon )$ is a $S_{j}$-symmetric periodic solution of the Hamiltonian system associated with the Hamiltonian (1.1) or (2.1) with fixed period $T=2\pi L_0^3$.

(ii)   Moreover, under the assumptions (a) and (b) in (2.12), there are $2$-parameter families of initial conditions

\[ {\bf Y}_{\delta L, \epsilon}^{j,k}={\bf Y}_{0}^{j,k}+{\bf Y}^{j,k}(\delta L,\epsilon)\quad {\bf Y}^{j,k}(\delta L,\epsilon)=(0,0,\delta L^{j,k},\delta G^{j,k}(\delta L,\epsilon)), \]

parametrized by $\epsilon$ and $\delta L$ sufficiently small, such that $\varphi (t,{\bf Y}_{\delta L, \epsilon }^{j,k};\epsilon )=\varphi _{0}(t, {\bf Y}_{0}^{j,k})+{\mathcal {O}}(\epsilon )$ is a $S_{j}$-symmetric periodic solution of the Hamiltonian system associated with the Hamiltonian (1.1) or (2.1) with differentiable period $\overline {T}(\delta L,\epsilon )=2\pi L_0^3+{\mathcal {O}}(\epsilon ^{\alpha })$.

Furthermore, if $\lambda _1, \lambda _2$ are the eigenvalues of $A$ in (2.11), then the characteristic multipliers of any of the periodic solutions $\varphi (t,{\bf Y}_{\epsilon }^{j,k};\epsilon )$ or $\varphi (t,{\bf Y}_{\delta L,\epsilon }^{j,k};\epsilon )$ are $1,1,1+\epsilon ^{\alpha } \lambda _1 +{\mathcal {O}}(\epsilon ^{\alpha +1}),1+\epsilon ^{\alpha }\lambda _2+{\mathcal {O}}(\epsilon ^{\alpha +1})$.

Proof. Our proof works for each case of symmetry, that is, $j=1$, or $j=2$. Let $\varphi _{0}(t, {\bf Y}^{j,k})$ be a solution of the Kepler problem as in (2.9) and $\varphi (t)= \varphi (t,{\bf Y}^{j,k};\epsilon )$ (with the same initial condition) be a solution of the perturbed Hamiltonian system associated with (2.1) in Delaunay variables. Since ${\bf Y}^{j,k} \in \mathcal {L}_{j}$, the solution $\varphi (t)$ will be $S_{j}$-symmetric, if at the instant $t=T/2$, it intercepts orthogonally the subspaces $\mathcal {L}_{j}$. By lemmas 2.1 and (2.4), we must verify the following two conditions:

\begin{align*} \ell(t,\delta L,\delta G,\epsilon) & = (L_0+\delta L)^{{-}3} \ t +{\mathcal{O}}(\epsilon^\alpha)= \pi,\\ g(t,\delta L, \delta G,\epsilon) & = g_{0}^{j,k}+\epsilon^{\alpha} g^{(1)}(t,{\bf Y}^{j,k}) +{\mathcal{O}}(\epsilon^{\alpha+1})= g_{0}^{j,k} +m\pi, \end{align*}

at the instant $t=T/2$ for some $m\in \mathbb {N}$. Taking $m=0$, the previous system (called periodicity equations) can be rewritten as

(2.13)\begin{equation} \begin{aligned} & f_1^j(t,\delta L,\delta G,\epsilon)= (L_0+\delta L)^{{-}3} \ t- \pi +{\mathcal{O}}(\epsilon^\alpha)= 0,\\ & f_2^j(t,\delta L, \delta G,\epsilon)= g^{(1)}(t,{\bf Y}^{j,k}) +{\mathcal{O}}(\epsilon)= 0, \end{aligned} \end{equation}

at the instant $t=T/2$. A natural but not trivial way to solve (2.13) is to apply the implicit function theorem. Under the choice of $T$ and hypothesis (a), it is clear that $f_1^j(T/2, 0, 0, 0)=f_2^j(T/2, 0, 0, 0)=0$ for $j=1,2$. Moreover, by differentiating the system (2.13) with respect to $(\delta L,\delta G)$ and evaluating at $t=T/2$, ${\bf Y}^{j,k}={\bf Y}_0^{j,k}$ and $\epsilon =0$, we obtain that the Jacobian matrix satisfies

\[ \frac{\partial (f_{1}^j,f_{2}^j)}{\partial (\delta L,\delta G)}\Big|_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_0^{j,k},\epsilon=0} = \left( \begin{array}{cc} -3 L_{0}^{{-}4}T/2 & 0 \\ \dfrac{\partial g^{(1)}}{\partial \delta L} & \dfrac{ \partial g^{(1)}}{\partial \delta G} \end{array}\right)_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0}. \]

It follows that

(2.14)\begin{equation} \det \frac{\partial (f_{1}^j,f_{2}^j)}{\partial (\delta L,\delta G)}\Big|_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0}={-} \frac{3 T}{2 L_0^4} \displaystyle\frac{ \partial g^{(1)}}{\partial \delta G}\Big|_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0}\neq0, \end{equation}

because by hypothesis (b)

\[ \frac{ \partial g^{(1)}}{\partial \delta G}\Big|_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0}= {\displaystyle \int_{0}^{T/2}}\frac{\partial}{\partial \delta G}\left(\frac{\partial {\mathcal{H}}_{1}}{\partial G }\right)(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf {\bf Y}^{j,k}}={\bf Y}_{0}^{j,k}}\,{\rm d}\tau\neq 0. \]

Thus, by the implicit function theorem, we obtain unique differentiable functions $\delta L^{j,k}=\delta L(\epsilon )$ and $\delta G=\delta G(\epsilon )$ defined for $\epsilon$ sufficiently small, such that $\delta L(0)=0$, $\delta G (0)=0$ and $f_i^j(T/2,\delta L(\epsilon ),$ $\delta G(\epsilon ),\epsilon )=0$ for $i,j=1, 2$. Thus, we obtain a $1$-parameter (on $\epsilon$) family of initial conditions

(2.15)\begin{equation} {\bf Y}_{\epsilon}^{j,k}=(0,g_0^{j,k},L_0+\delta L^{j,k}(\epsilon),G_0+\delta G^{j,k}(\epsilon)) \end{equation}

such that it gives rise to $S_j$-symmetric periodic solutions of the perturbed problem (2.1) with fixed period $T=2\pi L_0^3$. This proves item (i).

To prove item (ii), we introduce the time as a new independent variable in system (2.13). Again, for the hypothesis (a) it is clear that $f_1^j(T/2, 0, 0, 0)=f_2^j(T/2, 0, 0, 0)=0$. Moreover, by differentiating the system (2.13) with respect to $(t,\delta G)$ and evaluating at ${\bf Y}^{j,k}={\bf Y}_{0}^{j,k}$, $t=T/2$ and $\epsilon =0$, after some calculations, we verify that the Jacobian matrix satisfies

\[ \frac{\partial (f_{1}^j,f_{2}^j)}{\partial (t,\delta G)}\Big|_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0} = L_0^{{-}3} \frac{\partial g^{(1)}}{\partial \delta G}_{t=\frac{T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon=0}. \]

Again the hypothesis (b) implies that $\det \frac {\partial (f_{1}^j,f_{2}^j)}{\partial (t,\delta G)}\Big |_{t=\frac {T}{2},{\bf Y}^{j,k}={\bf Y}_{0}^{j,k},\epsilon =0}\neq 0$. Thus, by the implicit function theorem, we obtain unique differentiable functions $\delta G=\delta G(\delta L,\epsilon )$ and $\tau =\tau (\delta L,\epsilon )=T/2+ {\mathcal {O}}(\epsilon ^{\alpha })$ defined for $\epsilon$ and $\delta L$ sufficiently small, such that $\delta G(0,0)=0$ and $\tau (0,0)=T/2$ and $f_i^j(T/2,\delta L, \delta G(\delta L,\epsilon ),\epsilon )=0$ for $i, j=1, 2$. Therefore, we obtain a periodic $S_{j}$-symmetric solution of the perturbed system associated with the Hamiltonian function (2.1) with initial condition

(2.16)\begin{equation} {\bf Y}_{\delta L,\epsilon}^{j,k}=(0,g_0^{j,k},L_0+\delta L^{j,k},G_0+\delta G^{j,k}(\delta L^{j,k},\epsilon)), \end{equation}

which is $2\tau$-periodic and close to $T=2\pi L_0^3$ which is $\bar {T}=2\tau$-periodic such that $\bar {T}=2\pi L_0^3+\mathcal {O}(\epsilon ^{\alpha })$. Thus, we have proved item (b).

To study the type of stability of the previous periodic solutions, we are going to calculate its characteristic multipliers. Let

\[ \Sigma=\{(\ell,g,L,G) \ : \ {\mathcal{H}}(\ell,g,L,G)=h_0, \, \ell=0\}, \]

be a local cross-section on the level ${\mathcal {H}}=-\frac {1}{2L_0^2}= h_0$ in a neighbourhood of the point $(Y_1^{(0)}, Y_2^{(0)},$ $Y_3^{(0)}, Y_4^{(0)})$ given in (2.7). We denote by $X=(X_1,X_2)$ the points in $\Sigma$. Thus, considering $\overline {Y}= (Y_1^{(0)}, Y_2^{(0)}+ X_1, Y_3^{(0)}, Y_4^{(0)}+X_2)$ the Poincaré map $P$ on $\Sigma$ is given by $P(X,\epsilon )=(g(\mathcal {T},\overline {Y},\epsilon ),G(\mathcal {T},\overline {Y},\epsilon )),$ where $g$ and $G$ were characterized in (2.4) and $\mathcal {T}$ is the return time which is close to $T$. Using the form of ${\bf Y}_{0}^{j,k}$ and (2.5)–(2.6), we arrive at

(2.17)\begin{align} P(X,\epsilon)& = (g_{0}^{j,k}+X_1, G_0+X_2)\nonumber\\ & \quad +\epsilon^{\alpha} \left(\displaystyle\int_{0}^{\mathcal{T}}\frac{\partial {\mathcal{H}}_1}{\partial G}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t, -\displaystyle\int_{0}^{\mathcal{T}}\frac{\partial {\mathcal{H}}_1}{\partial g}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t\right)+ \mathcal{O}(\epsilon^{\alpha+1}),\nonumber\\ & = (g_{0}^{j,k}+X_1,G_0+ X_2)\nonumber\\ & \quad +\epsilon^{\alpha} \left(\int_{0}^{T}\frac{\partial {\mathcal{H}}_1}{\partial G}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t, -\int_{0}^{T}\frac{\partial {\mathcal{H}}_1}{\partial g}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t\right)+\mathcal{O}(\epsilon^{\alpha+1}). \end{align}

Recalling (2.10), we have that the differential of $P$ has the form

(2.18)\begin{equation} DP(X, \epsilon)= I+ \epsilon^{\alpha} D_X \overline{H}(X)+ \mathcal{O}(\epsilon^{\alpha+1}). \end{equation}

Since the $1$-parameter family of initial conditions of the $T$-symmetric periodic solutions are $Y_{0}^{j,k}+ (0, 0, \delta L(\epsilon ), \delta G(\epsilon ))$, then the respective points on the section cross $\Sigma$ will be $X_{\epsilon }= ( g_{0}^{j,k}, G_0+\delta G(\epsilon ))$ with $X_0=( g_{0}^{j,k}, G_0)$. Therefore,

(2.19)\begin{equation} DP(X_{\epsilon}, \epsilon)= I+ \epsilon^{\alpha} A+ \mathcal{O}(\epsilon^{\alpha+1}), \end{equation}

with $A$ as in (2.11). Since, the nontrivial characteristic multipliers of the $T$-symmetric periodic solutions are the eigenvalues of $DP(X_{\epsilon }, \epsilon )$, we obtain the result to the periodic solutions generated by the $1$-parameter family of initial conditions. To study the type of stability for periodic solutions given by the $2$-parameter families of initial conditions, we follow the same ideas as in the previous case. Again, the return map $P$ is given by $P(X, \epsilon )=(g(\mathcal {T}, \overline {Y},\epsilon ),G(\mathcal {T} ,\overline {Y},\epsilon ))$, where $\mathcal {T}$ is the return time which is close to $\overline {T}(\delta L, \epsilon )= 2\pi L_{0}^{3}+{\mathcal {O}}(\epsilon ^{\alpha })$. From (2.5) to (2.6), it follows that

(2.20)\begin{align} P(X,\epsilon)& = (g_{0}^{j,k}+X_1, G_0+X_2)\nonumber\\ & \quad +\epsilon^{\alpha}\left(\displaystyle\int_{0}^{\mathcal{T}}\frac{\partial { H}_1}{\partial G}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t, -\int_{0}^{\mathcal{T}}\frac{\partial {H}_1}{\partial g}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t\right)+ \mathcal{O}(\epsilon^{\alpha+1}),\nonumber\\ & = (g_{0}^{j,k}+X_1, G_0+X_2)\nonumber\\ & \quad +\epsilon^{\alpha} \left(\int_{0}^{\overline{T}}\frac{\partial { H}_1}{\partial G}(\varphi_{0}(t,\overline{Y}))\,{\rm d}t, -\displaystyle\int_{0}^{\overline{T}}\frac{\partial {H}_1}{\partial g}(\varphi_{0}(t,\overline{Y}))\,{\rm d}t\right)+ \mathcal{O}(\epsilon^{\alpha+1})\nonumber\\ & = (g_{0}^{j,k}+X_1, G_0+X_2)\nonumber\\ & \quad +\epsilon^{\alpha} \left(\displaystyle\int_{0}^{T}\frac{\partial { H}_1}{\partial G}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t, -\displaystyle\int_{0}^{T}\frac{\partial {H}_1}{\partial g}(\varphi_{0}(t, \overline{Y}))\,{\rm d}t\right)+ \mathcal{O}(\epsilon^{\alpha+1}). \end{align}

Since the initial condition of the $\bar {T}$-symmetric periodic solutions is $Y_{0}^{j,k}+$ $(0,0, \delta L, \delta G(\delta L,\epsilon ))$, then the respective points on the local cross-section $\Sigma$ will be $X_{\epsilon }= (g_{0}^{j,k}, G_0+\delta G(\delta L,\epsilon )))$ with $X_0=(g_{0}^{j,k},G_0)$. Then, the nontrivial characteristic multipliers associated with the symmetric $\overline {T}$-periodic solutions $\varphi (t,{\bf Y}_{\delta L, \epsilon }^{j,k};\epsilon )$, given by the $2$-parameter families of initial conditions, are the eigenvalues of (2.19). Thus, we have proved the theorem.

Remark 2.4 It is clear that we cannot obtain double-symmetric periodic solutions as continuation of elliptic Keplerian solutions, unless it is circular.

Remark 2.5 Note that fixed $j$ and the energy level $\mathcal {H}_0= -\frac {1}{2L_0^2}<0$ of the Keplerian solution, then for each $G_0$ (in fact, $G_0(L_0)$, i.e., $G_0$ is in general a function of $L_0$) solution of equation (a) of theorem 2.3 (such that (b) is verified) gives us two families (distinct) of initial conditions such that each of them arises $S_j$-symmetric periodic solution, one for $k=1$ and the other for $k=2$. In fact, introducing the invariants $a=(a_1, a_2, a_3)$ space according to [22], it is known that different points $a$ give us different orbits of the full system. For this purpose, let ${\bf A}$ be the Laplace–Runge–Lenz vector, $L_0=(-2h^*)^{-1/2}$ and $a={\bf G}+L_0{\bf A}$. One can check that $\|a\|=L_0$ and the vector $a$ determines uniquely an orbit of the Kepler problem on the energy level $h^*$. Each point of $\|a\|=L_0$ with $\|a\|\neq G_0\neq 0$ corresponds to an elliptic orbit of the Kepler problem. Since, explicitly $a_1=eL_0\sin g$, $a_2=eL_0\cos g$ and $a_3=G_0$ evaluating them in $g=g_0^{jk}$, we obtain different vectors $a$ for $k=1$ and $k=2$ and $j$-fixed. Thus, we have proved the affirmation.

The main difficulty imposed by the conditions (a) and (b) in theorem 2.3 is related to the calculus of the partial derivatives and the integration. In some situations, for practical problems, the main perturbation term $H_1$ cannot be obtained in a closed-form in Delaunay variables. One strategy in order to compute the partial derivatives involved in the conditions given by theorem 2.3 consists of writing the perturbed function $H_1$ in mixed coordinates which involves the polar and Delaunay elements (see [22] and references therein for more details about these coordinates). For the integration, we can make an appropriate change of variables using the Kepler equation so we integrate with respect to the eccentric anomaly. More precisely, to compute the partial derivative $\frac {\partial \mathcal {H}_1}{\partial G}$, we consider the Kepler equation

(2.21)\begin{equation} \ell=E-e\sin E, \end{equation}

where $e=\sqrt {1-G^2/L^2}$ is the eccentricity of the Keplerian orbit and $E$ is the eccentric anomaly.

The perturbed function $\mathcal {H}_1$, in the mixed polar and Delaunay elements, depends on the variables $r$, $f$, $g$, $L$ and $G$. First, we eliminate the dependence of $\mathcal {H}_1$ on the variables $r$ and $f$, using the auxiliary relations

(2.22)\begin{equation} r=\frac{a(1-e^2)}{1+e \cos f}, \quad \cos f=\frac{\cos E -e}{1-e \cos E}, \quad \sin f=\frac{\sqrt{1-e^2}\sin E}{1-e \cos E}\quad e=\sqrt{1-G^2/L^2}, \end{equation}

where $a=L^2$. After these substitutions, the function $\mathcal {H}_1$ is simply a function of the form

(2.23)\begin{equation} \mathcal{H}_1=\mathcal{H}_1\left(E(G),g,L,G\right). \end{equation}

Since $\frac {\partial e}{\partial G}=-G/(L^2e)$, by differentiation of (2.21) with respect to $G$, we obtain

(2.24)\begin{equation} \frac{\partial E}{\partial G}={-}\frac{G\sin E}{L^2 e(1-e\cos E)}. \end{equation}

Thus, $\frac {\partial \mathcal {H}_1}{\partial G}$ can be calculated using equation (2.23) and the chain rule. To integrate equations (2.12), we introduce the change in the variables

(2.25)\begin{equation} t=L_0^3(E-e_0 \sin E), \end{equation}

where $L_0^2$ and $e_0$ correspond to the semi-major axis and the eccentricity of the Keplerian orbit $\varphi _{0}(\tau, {\bf Y}_0^{j,k}))$, respectively.

2.1 Relationship with the averaging method

Here we intend to find some relations or comparisons between the periodic solutions obtained by the continuation of the elliptic solutions (non-circular) of the Kepler problem given in theorem 2.3 with those found by using the averaging theory on Hamiltonian systems.

Initially, in order to apply the averaging method, let $h^*$ be a negative real constant such that $\mathcal {H}_0(L)=-\frac {1}{2L^2}=h^*$, so $\mathcal {H}_0^{-1}$ in a neighbourhood of $h^*$ is a diffeomorphism and $\mathcal {H}_0^{-1}(h^*)$ is a compact connected circle bundle over the base space $B(h^*)$ with projection $\pi :\mathcal {H}_0^{-1}(h^*)\rightarrow B(h^*)$ (see [36] for more details). Moreover, all the solutions of the Hamiltonian system associated with the Hamiltonian (2.2) are periodic and have periods depending smoothly on $h^*$, i.e., the period is a smooth function $T=T(h^*)$.

Next, we denote by $\overline {\mathcal {H}}$ the averaging function with respect to the mean anomaly $\ell$ as

(2.26)\begin{equation} \overline{\mathcal{H}}=\frac{1}{2\pi}\int_0^{2\pi}\mathcal{H}_1(\ell,g,L,G)d\ell=\frac{1}{2\pi}\int_0^{2\pi}\mathcal{H}_1(E-e\sin E,g,L,G)L(1-e\cos E)\,{\rm d}E. \end{equation}

Thus, at the energy level $\mathcal {H}=h^*$ ($h^*<0$), the averaged system associated with (2.26) is given by

(2.27)\begin{equation} \begin{aligned} \frac{d g}{d \ell} & =\epsilon \frac{1}{({-}2h^*)^{3/2}}\frac{\partial\overline{\mathcal{H}} }{\partial G}=\epsilon \frac{1}{({-}2h^*)^{3/2}}F_1(g, G),\\ \frac{d G}{d \ell} & ={-}\epsilon \frac{1}{({-}2h^*)^{3/2}}\frac{\partial\overline{\mathcal{H}} }{\partial g}=\epsilon \frac{1}{({-}2h^*)^{3/2}}F_2(g, G). \end{aligned} \end{equation}

Next, for our problem, we have the following result as a consequence of averaging theory or Reeb's theorem for Hamiltonian systems [22, 30].

Theorem 2.6 Consider the averaged differential system (2.27) in Delaunay variables restricted to the energy level $\mathcal {H}= h$ ($h<0$). If $\bar {p}=(g_0,G_0)$ is a non-degenerate critical point, then there are smooth functions $p(\epsilon )= (\ell _0(\epsilon ), g(\epsilon ), L(\epsilon ), G(\epsilon ))$ and $T(\epsilon )$ for $\epsilon$ small with $p(0)=(\ell _0,g_0,1/\sqrt {-2h^*},G_0)$ and $T(0)=T=2\pi L_0^3$ such that the solution of (2.1) through $p(\epsilon )$ is $T(\epsilon )$-periodic.

In addition, if the characteristic exponents of the critical point $\bar {p}$ (i.e., the eigenvalues of the matrix $A=\mathbb {J}D^2\mathcal {\bar {H}}(\bar {p})$) are $\lambda _1,\lambda _2$, then the characteristic multipliers of the periodic solution through $p(\epsilon )$ are

\[ 1,1,1+\epsilon\lambda_1 T+O(\epsilon^2),1+\epsilon\lambda_2 T+O(\epsilon^2). \]

Next, we will discuss the relationship between both methods, theorem 2.3 using Poincaré continuation and theorem 2.6 using the averaging theory. Our explanation is performing under two points of view: the curve of initial conditions and the properties of the family of periodic solutions obtained by each theorem.

Without loss of generality, we will assume that the Hamiltonian (1.1) [resp. Hamiltonian (2.1)] is invariant with respect to the symmetry $S_2$, i.e., the reflection with respect to the $x$-axis. The analysis of the other symmetry $S_1$ follows using similar arguments.

Since $\ell =E-e\sin E$, then the function $\mathcal {H}_1(\ell,g,L,G)$ assumes the form $\mathcal {H}_1=\mathcal {H}_1 (E,g,L,G)$. By the invariance of the function $\mathcal {H}_1$ with respect to the symmetry $S_2$ it follows that

(2.28)\begin{equation} \mathcal{H}_1(E,g,L,G)=\mathcal{H}_1(2\pi-E,2\pi-g,L,G).\end{equation}

Differentiating both sides of the last equation with respect to $g$ and evaluating at $g=k\pi$, $k=0,1$, we obtain

(2.29)\begin{equation} \frac{\partial \mathcal{H}_1}{\partial g}(E,k\pi,L,G)={-}\frac{\partial \mathcal{H}_1 }{\partial g}(2\pi-E,k\pi,L,G) \quad k=0,1. \end{equation}

Multiplying both sides of equation (2.29) by $L^3(1-e\cos E)$ and integrating in $E$ from $0$ to $2\pi$, we obtain

(2.30)\begin{align} & \int_0^{2\pi}\frac{\partial \mathcal{H}_1 }{\partial g}(E,k\pi,L,G)L^3(1-e\cos E)\,{\rm d}E\nonumber\\ & \quad ={-}\int_0^{2\pi}\frac{\partial \mathcal{H}_1 }{\partial g}(2\pi-E,k\pi,L,G)L^3(1-e\cos E)\,{\rm d}E. \end{align}

Making the change $\bar {E }=2\pi -E$ on the right-hand side of equation (2.30), we obtain (after removing the bar)

\[ \int_0^{2\pi}\frac{\partial {\mathcal{H}}_1}{\partial g}(E,k\pi,L,G)L^3(1-e\cos E)\,{\rm d}E=0. \]

Now, from (2.26), we have

(2.31)\begin{equation} \frac{\partial\overline{\mathcal{H}} }{\partial g}|_{g=k\pi}=\int_0^{2\pi}\frac{\partial {\mathcal{H}}_1}{\partial g}(E,k\pi,L,G)L^3(1-e\cos E)\,{\rm d}E=0, \quad k=0,1. \end{equation}

Thus, $F_2(g=k\pi, G)=0$ for $k=0,1$. In addition, from equation (2.31), it follows that $\frac {\partial ^2 \overline {\mathcal {H}}}{\partial G\partial g}\Big |_{g=k\pi }=0$.

Proposition 2.7 Suppose that the Hamiltonian function $\mathcal {H}$ given in (2.1) is invariant with respect to the symmetry $S_2$. Thus,

(2.32)\begin{equation} g^{(1)}(t,{\bf Y}_0^{2,k})=L_0^3\pi F_1(g=k\pi,G_0). \end{equation}

Proof. For $\mathcal {H}_1$ given in (2.23) we have

(2.33)\begin{align} g^{(1)}(t,{\bf Y}^{2,k})\Big|_{g=k\pi}& = \int_0^{\pi}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(E,g,L,G)\Big|_{g=k\pi}L^3(1-e\cos E)\,{\rm d}E\nonumber\\ & =L^3\frac{\partial }{\partial G}\left(\displaystyle\int_0^{\pi}{\mathcal{H}}_{1}(E,g,L,G)(1-e\cos E)\Big|_{g=k\pi}\,{\rm d}E\right)\nonumber\\ & =\frac{L^3}{2}\frac{\partial }{\partial G}\left(\int_0^{2\pi}{\mathcal{H}}_{1}(E,g,L,G)(1-e\cos E)\Big|_{g=k\pi}\,{\rm d}E\right)\nonumber\\ & =\frac{L^3}{2}\left(\frac{\partial }{\partial G}\int_0^{2\pi}{\mathcal{H}}_{1}(E,g,L,G)(1-e\cos E)\,{\rm d}E\right)\Big|_{g=k\pi}\nonumber\\ & ={\pi L^3}\frac{\partial \overline{H}_1 }{\partial G}\Big|_{g=k\pi}={\pi L^3}F_1(g=k\pi,G) \end{align}

Note that the equality

\begin{align*} & 2\int_0^{\pi}{\mathcal{H}}_{1}(e,E,g,L,G)(1-e\cos E)\Big|_{g=k\pi}\,{\rm d}E\nonumber\\ & \quad =\int_0^{2\pi}{\mathcal{H}}_{1}(e,E,g,L,G)(1-e\cos E)\Big|_{g=k\pi}\,{\rm d}E, \end{align*}

follows from the fact that by equation (2.28), the function $F(E)={\mathcal {H}}_{1}(e,E,g,L,G)(1-e\cos E)\Big |_{g=k\pi }$ satisfies $F(E)=F(2\pi -E)$. Finally, from (2.33) we obtain (2.32).

From the previous discussion, we can relate the $S_2$-symmetric periodic solutions obtained by theorem 2.3 with those obtained using the averaging method described in theorem 2.6. More precisely, from theorems 2.3 and 2.6 we obtain the following result.

Theorem 2.8 Consider $T=2\pi L_0^3$ and let $\varphi _{0}(t, {\bf Y}_{0}^{j,k})=(\frac {t}{L_{0}^3}, g_{0}^{j,k}, L_0, G_0)$ be a $S_j$-symmetric elliptic Keplerian solution $T$-periodic such that the conditions (a) and (b) of theorem 2.3 are satisfied.

Thus, $\bar {p}=(g_{0}^{j,k},G_0)$ is a critical point of the differential system (2.27). Furthermore, if $\bar {p}$ is an isolated critical point of (2.27) and $\frac {\partial ^2 \overline {\mathcal {H}}}{\partial g^2}(\bar {p})\neq 0$, then the elliptic Keplerian orbit $\varphi _{0}(t,{\bf Y}_0^{j,k})$ $($associated with the critical point $\bar {p})$ can be continued by Reeb's theorem to a $p(t, \bar {p}(\epsilon ))$, $T(\epsilon )$-periodic solution.

Moreover, if $1+\epsilon ^{\alpha } \lambda _1 +{\mathcal {O}}(\epsilon ^{\alpha +1}),1+\epsilon ^{\alpha }\lambda _2$ are the nontrivial characteristic multipliers of the periodic solutions $\varphi (t,{\bf Y}_{\epsilon }^{j,k};\epsilon )$ given by theorem 2.3, then the nontrivial multipliers characteristic of the periodic solutions $p(t,\bar {p}(\epsilon ))$ are

\[ 1+\epsilon^{\alpha} \lambda_1T +{\mathcal{O}}(\epsilon^{\alpha+1}),1+\epsilon^{\alpha}\lambda_2T. \]

Proof. We give the proof for the $S_2$-symmetric periodic solutions. Let $\varphi _{0}(t,{\bf Y}^{2,k})$ be an elliptic Keplerian solution near the elliptic solution $\varphi _{0}(t,{\bf Y}_0^{2,k})=(\frac {t}{L_{0}^3},0, L_0, G_0)$ such that the conditions (a) and (b) of theorem 2.3 are satisfied. By the previous discussion and by proposition 2.7 it follows that $\bar {p}=(g_0=k\pi,G_0)$ is a critical point of the average function $\bar {\mathcal {H}}$ given in (2.26).

Since $\frac {\partial ^2 \bar {\mathcal {H}}}{\partial g\partial G}\Big |_{g=k\pi }=0$, the non-degeneracy condition at the critical point $\bar {p}$ is given by

\[ \mathcal{D}=\det \mathbb{J}D^2\mathcal{H}(\bar{p})=\frac{\partial^2 \bar{\mathcal{H}}}{\partial G^2}\frac{\partial^2 \bar{\mathcal{H}}}{\partial g^2}\Big|_{g=k\pi,G=G_0}\neq 0. \]

From (2.32) it follows that

\[ \frac{\partial^2 \bar{\mathcal{H}}}{\partial G^2}(\bar{p})=\frac{1}{\pi L_0^3}\frac{\partial g^{(1)}(t,{\bf Y}_0^{2,k})}{\partial G}(\bar{p}), \]

and by condition (b) of theorem 2.3, we obtain $\frac {\partial ^2 \bar {\mathcal {H}}}{\partial G^2}(\bar {p})\neq 0$. Finally, since $\frac {\partial ^2 \overline {\mathcal {H}}}{\partial g^2}(p)\neq 0$ we obtain $\mathcal {D}\neq 0$. The conclusion about the nontrivial multipliers characteristic is immediate from theorems 2.3 and 2.6. Thus, we have concluded the proof.

We finish this section with some remarks about the similarities/differences of the periodic solutions obtained in theorem 2.3 and those given in theorem 2.8.

• • Let $\varphi _{0}(t,{\bf Y}_0^{2,k})=\left (\frac {t}{L_{0}^3},0, L_0, G_0\right )$ be an elliptic Keplerian solution and assume that the conditions (a) and (b) of theorem 2.3 are satisfied. By the previous discussion it follows that $\bar {p}=(g_0=0,G_0)$ is not necessarily a non-degenerate critical point of the system (2.26). Thus, the condition (b) of theorem 2.3 is a weaker condition when compared with the condition for the non-degeneracy of the critical point $\bar {p}$ given by theorem 2.8. For a concrete example of this situation, see remark 3.5.

• • Note that the periodic ($S_2$-symmetric) solutions obtained by theorem 2.3 item (i) have the same period $T=2\pi L_0^{3}$ that the elliptic orbit $\varphi _{0}(t,{\bf Y}_0^{2,k})$ of the unperturbed system. The solutions obtained in theorem 2.8 have period $T=2\pi L_0^{3}+\mathcal {O}(\epsilon )$.

• • The approximation of the period of the periodic solutions given in theorems 2.8 and 2.3 item (ii) can be computed as follows. If $\ell (0)=0$ and $\ell (T)=2\pi$, then $T(\varepsilon )$ must satisfy the equation

  \begin{align*} 2\pi & = \ell(T)-\ell(0)=\displaystyle\int_0^T \dot{\ell}(t)\,{\rm d}t= L_0^{{-}3} T+\varepsilon^{\alpha} L_0^{{-}3} \displaystyle\int\limits_0^{2 \pi} \dfrac{\partial \mathcal{H}_1} {\partial L}(g,L,G)\,{\rm d}\ell +O(\varepsilon^{\alpha+1})\\ & = L_0^{{-}3} T+ 2 \pi L_0^{{-}3} \varepsilon^{\alpha} \dfrac{\partial \mathcal{H}_1} {\partial L}(g,L,G)|_{L=L_0, g=g_0, G=G_0} +O(\varepsilon^{\alpha+1})\\ & = L_0^{{-}3} \left(T+ 2 \pi\varepsilon^{\alpha} \dfrac{\partial \bar{\mathcal{H}}} {\partial L}(g,L,G)|_{L=L_0, g=g_0, G=G_0}\right) +O(\varepsilon^{\alpha+1}) . \end{align*}

  By use of the implicit function theorem, it follows that $T(\varepsilon )=2\pi L_0^3+\varepsilon T^\ast +O(\varepsilon ^2)$, where $T^\ast =-2 \pi \dfrac {\partial \bar {\mathcal {H}}}{\partial L}(g,L,G)|_{L=L_0, g=g_0, G=G_0}$.

  We call attention for the initial conditions generating periodic solutions of the full system on the reduced space $\mathcal {B}(h)$ (parametrized by the coordinates $(g, G)$) obtained in theorem 2.3 and those given by theorem 2.8. The curve of initial conditions parametrized by $\epsilon$ at a fixed point $\bar {p}=(g_0^{j,k},G_0)$ as in corollary 2.8 has the form $(g(\epsilon ), G(\epsilon ))$. On the other hand, the curve obtained in theorem 2.3 at the same point $\bar {p}$ has the form $(k \pi, \overline {G}(\epsilon ))$. In figure 3, we illustrate these two families of initial conditions in a neighbourhood of a fixed critical point $\bar {p}$ on the reduced space.

• • The periodic solutions given by theorem 2.8 are not necessarily symmetric. The authors in reference [36] (see § 2.4) did the inverse process of this section, i.e., they showed that the Keplerian solutions associated with the critical point of the averaging system can be continued to a symmetric periodic solution under some restrictions. In fact, the technique of these authors consists in the use of the symmetries of the problem combined with averaging, reduction theory and the use of the implicit function theorem of Arenstorf. But, the use of the theorem of Arenstorf requires a suitable choice of the parameter $\epsilon$. In fact, the value of the parameter $\epsilon$ must be chosen in a discrete set. Thus, the symmetric periodic solutions analysed in [36] have the following properties: the existence of the symmetric solutions is conditioned only for some values of the parameter $\epsilon$; the period of these symmetric periodic solutions is very large and is given by $\tau (\epsilon )=\beta T+\mathcal {O}(\epsilon )$ (here, the parameter $\beta$ is a positive and very large integer and $T$ is the period of the Keplerian orbit).

  

  Figure 3. Illustration of the projection curves of initial conditions in a neighbourhood of the point $\bar {p}=(g_0=\pi, G_0)$ parametrized by the coordinates $(g, G)$. The blue curve $(g(\epsilon ), G(\epsilon ))$ represents the curve of initial conditions obtained by theorem 2.8. The red curve $( \pi, \overline {G}(\epsilon ))$ is the projection of the curve of initial conditions given by theorem 2.3 over the reduced space and it is associated with the $S_2$-symmetric periodic solutions.

  The last paragraph shows that there are important differences between the symmetric solutions obtained in [36] and those obtained in § 2.4, since these can be obtained for all values of $\epsilon$ (sufficiently small) and there is no restriction on the period $T$.

3. Applications

3.1 Planar hydrogen atom with Stark and quadratic Zeeman effect

We consider a simple atomic system, namely hydrogen atom interacting with time-independent and external fields that include electric and magnetic fields. In this formulation, the magnetic field is usually applied perpendicular to the $xy$-plane, while the static electric field is applied along the $x$-axis. See [11, 12, 15] or [24] for more details on the formulation of this problem. We call this problem simply by HCEM problem, whose Hamiltonian has the form

(3.1)\begin{equation} H(x,y,p_x,p_y)=\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{1}{\sqrt{x^2+y^2}}+\frac{\tilde{B}}{2}(yp_x-xp_y)+\frac{B^2}{8}(x^2+y^2)-F x. \end{equation}

This Hamiltonian function depends on two parameters $\tilde {B}$ and $F$. The term $Fx$ is the electrostatic potential describing the Stark effect while the other terms having the parameter $\tilde {B}$ refer to the linear and quadratic Zeeman effect (see [12, 24] for more details for definition of the physical constants). An important point in this problem is that the angular momentum is not a conserved quantity. We pint out that the Hamiltonian function (3.1) is invariant under (only) the symplectic reflection $S_2:(x,y,p_x,p_y)\rightarrow (x,-y,-p_x,p_y)$. We will consider the hydrogen atom in crossed electric and magnetic field problem for weak magnetic and electric fields. This is achieved scaling $\tilde {B}$ and $F$ in the following way:

\[ \tilde{B}=2 B \epsilon, \quad F=\gamma \epsilon, \]

where $\epsilon$ is a parameter sufficiently small. Considering this scaling, the Hamiltonian function (3.1) assumes the form

(3.2)\begin{equation} H(x,y,p_x,p_y)=\displaystyle\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{1}{\sqrt{x^2+y^2}}+\epsilon\left[B(yp_x-xp_y)-\gamma x\right]+\epsilon^2\frac{B}{2}(x^2+y^2). \end{equation}

To show the existence of symmetric periodic solutions, first as in the previous case, we write the Hamiltonian function (3.2) in mixed variables involving polar and Delaunay variables, so obtaining the Hamiltonian function

\[ \mathcal{H}={-}\frac{1}{2L^2}+\epsilon\mathcal{H}_1(E,g,L,G)+\mathcal{O}(\epsilon ^2), \]

where

(3.3)\begin{align} \mathcal{H}_1(E,g,L,G)& ={-}B G - r \gamma \cos(f+g),\nonumber\\ & ={-}B G+a \gamma \left(\sqrt{1-e^2} \sin E \sin (g)+\cos (g) (e-\cos E)\right), \end{align}

$r=a(1-e^2)/(1+e\cos f)$ and $e=\sqrt {1-(G^2/L^2)}$. The main result about the existence of $S_2$-symmetric periodic solutions for the $HCEM$ problem (3.2) is the following.

Theorem 3.1 Fix the energy level ${\mathcal {H}}_{0}=- \frac {1}{2 L_0^2}$ and the period $T=2\pi L_0^3$ of the elliptic Kepler solution. Consider $\gamma$ and $B_1=-B/\gamma$ non null real constants and $\epsilon$ sufficiently small. Thus, for the planar HCEM problem (3.2), there exist two 1-parameter (and two 2-parameter) families of initial conditions such that each of them gives us a second-kind $S_2$-symmetric periodic solutions. Two of these initial conditions give us prograde solutions and the other two retrograde solution. All the periodic symmetric solutions are linearly stable.

Proof. Maintaining the notation of § 2, first we use expression (2.24) for calculations involving the derivative $\frac {\partial {\mathcal {H}}_{1}}{\partial G}$. After some calculations, we obtain

(3.4)\begin{equation} \frac{\partial {\mathcal{H}}_{1}}{\partial G}={-}B+G \gamma \cos g \left(-\frac{\sin ^2E}{e-\text{ee}^2 \cos E}-\frac{1}{e}\right)-\frac{G \gamma \sin E \sin g (e-\cos E)}{e \sqrt{1-e^2} (e \cos E-1)}. \end{equation}

Next, we evaluate it in the solution $\varphi _{0}(\tau, {\bf Y}^{2,k})$, and after some simplifications, we arrive at

(3.5)\begin{align} & \frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{2,k})) \nonumber\\ & \quad =\frac{B e+({-}1)^kG \gamma-e \cos (E) (B e+({-}1)^{k+1}G \gamma )+({-}1)^kG \gamma \sin ^2(E) }{e (e \cos (E)-1)},\end{align}

where $e=\sqrt {1-G^2/L^2}$, $G=G_0+\delta G$, $L=L_0+\delta L$. For a fixed value of $0< e<1$ the function given by (3.5) is continuous and differentiable with respect to $E$ and then the integral of (3.5) can be calculated explicitly. For the integration of equation (3.5), we use the change of variables (2.25). After the integration of equation (3.5) and some simplifications, we arrive at

(3.6)\begin{align} g^{(1)}(T/2,{\bf Y}_0^{2,k})& = \displaystyle\int_{0}^{\pi}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{2,k}))|_{{\bf Y}^{2,k}={\bf Y}_{0}^{2,k}}L_{0}^{3}(1-e_0\cos E)\,{\rm d}E\nonumber\\ & =\pi L_0^3 \left(- B+({-}1)^{k+1}\frac{3 \gamma G_0 L_0}{2\sqrt{L_0^2-G_0^2}}\right)\nonumber\\ & =\pi L_0^3 \gamma \left(- \frac{B}{\gamma}+({-}1)^{k+1}\frac{3 G_0 L_0}{2\sqrt{L_0^2-G_0^2}}\right). \end{align}

To find the solutions (in the variable $G_0$) of $g^{(1)}(T/2,{\bf Y}_0^{2,k})=0$, we define the normalized parameter $B_1=-B/\gamma$. Moreover, we introduce the auxiliary functions

(3.7)\begin{equation} f_1^k(G_0)=({-}1)^{k+1}\frac{3 G_0 L_0}{2\sqrt{L_0^2-G_0^2}},\quad f_2(G_0)=B_1,\quad k=0,1, \end{equation}

defined in $(-L_0,L_0)$. Of course, $g^{(1)}(T/2,{\bf Y}_0^{2,k})=0$, if and only if,

(3.8)\begin{equation} f_1^k(G_0)={-}f_2(G_0), \end{equation}

for some $G_0\in (-L_0,L_0)$. Due to straightforward properties of the function $f_1^k$, we obtain the following possibilities for the solutions of (3.8):

• • If $k=0$ and $B_1<0$ (resp. $B_1>0$), then there is a unique solution $G_0^{(0)}=\frac {2 B_1 L_0}{\sqrt {4 B_1^2+9 L_0^2}}<0$ (resp. $G_0^{(0)}=\frac {2 B_1 L_0}{\sqrt {4 B_1^2+9 L_0^2}}>0$).

• • If $k=1$ and $B_1<0$ (resp. $B_1>0$), then there is a unique solution $G_0^{(1)}=-\frac {2 B_1 L_0}{\sqrt {4 B_1^2+9 L_0^2}}>0$ (resp. $G_0^{(1)}=-\frac {2 B_1 L_0}{\sqrt {4 B_1^2+9 L_0^2}}<0$).

Observe that $0<|G_0^{(0)}|=|G_0^{(1)}|< L_0$, whenever $B_1\neq 0$ and $\gamma \neq 0$. Thus, we have verified the condition (a) of theorem 2.3.

Now, we are going to verify the condition (b) given by theorem 2.3. Consider the equation (2.24) and the chain rule to differentiate (3.5) with respect to $\delta G$. After that, we evaluate $\displaystyle \frac {\partial }{\partial \delta G}\left (\frac {\partial \mathcal {H}_1}{\partial G}\right )$ on the solution $\varphi _{0}(\tau, {\bf Y}_0^{2,k}))$. Next, using the change (2.25), after integration and simplification, we arrive at

(3.9)\begin{align} \frac{\partial g^{(1)}}{\partial \delta G}\Big |_{{\bf Y}^{2,k}={\bf Y}_{0}^{2,k}} & =\displaystyle\int_{0}^{\pi}\frac{\partial }{\partial \delta G}\left(\frac{\partial \mathcal{H}_1}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{2,k}))\right)\Big |_{{\bf Y}^{2,k}= {\bf Y}_{0}^{2,k}}L_0^3(1-e_0\cos E)\,{\rm d}E \nonumber\\ & ={-}\frac{3 \pi \gamma L_0^6}{2 \left(L_0^2-G_0^2\right)^{3/2}}. \end{align}

Of course, for all $L_0,\gamma \neq 0$ we have $\frac {\partial g^{(1)}}{\partial \delta G}\Big |_{{\bf Y}_{0}^{2,k}}\neq 0$ and the item (b) of theorem 2.3 is verified. Thus, we obtain the existence of two $1$-parameter (on $\epsilon$) families of initial conditions

\[ {\bf Y}_{\epsilon}^{2,k}=\left(0, k\pi, L_0+\delta L^{(k)}(\epsilon), G_0^{(k)}+\delta G^{(k)}(\epsilon)\right), \quad k=0,1, \]

where $G_0^{(k)}=(-1)^{k}\frac {2 B_1 L_0}{\sqrt {4 B_1^2+9 L_0^2}}$, such that each of them gives us a $S_2$-symmetric of second-kind periodic solution with period $T=2\pi L_0^3$. In addition, we have two $2$-parameters (on $\epsilon$ and $\delta L$) families of initial conditions

\[ {\bf Y}_{\epsilon,\delta L}^{2,k}=(0, k\pi, L_0+\delta L^{(k)},G_0^{(k)}+\delta G^{(k)}(\epsilon,\delta L^{(k)})), \quad k=0,1, \]

such that, each of them gives us $S_2$-symmetric periodic solution of second kind with period $\overline {T}$ close to $T=2\pi L_0^3$.

Now, we discuss the stability of the previous periodic solutions. After some algebraic manipulation, the averaging function $\overline {H}$ given in (2.10) here has the form

(3.10)\begin{equation} \overline{H}(X_1,X_2)=\pi L_0^3 \left(3 \gamma L_0 \cos (g_0^{2,k}+X_1) \sqrt{L_0^2-(G_0+X_2)^2}+2 B (G_0+X_2)\right). \end{equation}

By the definition of the matrix $A$ in (3.10), we obtain

(3.11)\begin{equation} A=\left(\begin{array}{cc} 0 & -({-}1)^k\dfrac{3 \pi \gamma L_0^6}{\left(L_0^2-G_0^2\right)^{3/2}} \\ ({-}1)^k 3 \pi \gamma L_0^4 \sqrt{L_0^2-G_0^2} & 0 \\ \end{array}\right). \end{equation}

Thus, the corresponding characteristic multipliers of the prograde and retrograde $S_2$-symmetric periodic solutions are

\[ 1,1, 1+\epsilon \frac{3 i \pi \gamma L_0^5}{\sqrt{L_0^2-G_0^2}}+\mathcal{O}(\epsilon^2),1-\epsilon \frac{3 i \pi \gamma L_0^5}{\sqrt{L_0^2-G_0^2}}+ \mathcal{O}(\epsilon^2), \]

where $G_0=G_0^{(k)}$. Thus, in any case the $S_2$-symmetric periodic solutions are linearly stable.

3.2 Anisotropic two-body problem under Seeliger's potential

The gravitational potential of the two-body problem due to Seeliger's theory is given by

(3.12)\begin{equation} V({\bf q})={-}\frac{A}{\|{\bf q}\|}\,{\rm e}^{{-}K\|{\bf q}\|}, \end{equation}

where ${\bf q}$ is the vector between the two masses $m_1$ and $m_2$, $A$ and $K$ are positive constants. See [23, 26, 29] for more details on the formulation of this problem. We mention that a field featured by a (3.12)-like potential has larger physical implications. London's theory of superconductivity involves an electromagnetic potential of this form. Debye–Huckel's theory of screening in electrolytes leads to a similar screened potential. See [26] and references therein for more physical implications of this problem.

In this paper, following the formulation given in [26], we consider the Seeliger potential in an anisotropic space, where the potential function (3.12) depends on the parameter $\mu$ that measures the strength of anisotropy. More precisely, we consider the Hamiltonian

(3.13)\begin{equation} H(x,y,p_x,p_y)=\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{A}{\sqrt{\mu x^2+y^2}}e^{{-}K\sqrt{\mu x^2+y^2}}. \end{equation}

We call this problem as anisotropic Seeliger's Hamiltonian or shortly ASH. It is clear that the Hamiltonian function (3.13) is invariant under the symmetries $S_1$ and $S_2$.

To obtain a convenient approach of the Hamiltonian (3.13), we introduce the scaling $\mu =1-\mu _0\epsilon ^2$ and $K=\epsilon \xi$ for $\epsilon$ small. In addition, we make the change ${\bf {p}}=\sqrt [3]{A}{\bf {p}}$, ${\bf {q}}=\sqrt [3]{A}{\bf {q}}$, that is, a $1/A^{2/3}-$symplectic change. After this scaling, we develop the resulting Hamiltonian in a Taylor series in $\epsilon$ around $\epsilon =0$. Thus, eliminating the constant terms and setting $\kappa =\xi A^{1/3}$, ending up with

(3.14)\begin{align} H(x,y,p_x,p_y)& =\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{1}{\sqrt{x^2+y^2}}\nonumber\\ & \quad - \epsilon ^2 \left(\frac{1}{2} \kappa^2 \sqrt{x^2+y^2}+\frac{x^2\mu_0}{2 \left(x^2+y^2\right)^{3/2}}\right)+\mathcal{O}(\epsilon^3). \end{align}

Remark 3.2 Note that the Hamiltonian function (3.14) depends on two parameters $\kappa$ and $\mu _0$. The parameter $\mu _0$ determines the direction of the predominant force. In the case in which $\mu _0>0$ ($\mu <1$), the attraction is weakest in the direction of the $x$-axis and strongest in that of the $y$-axis. The situation is reversed if $\mu _0<0$ ($\mu >1$).

In the mixed polar and Delaunay coordinates, the Hamiltonian (3.14) assumes the form

(3.15)\begin{equation} \mathcal{H}={-}\frac{1}{2L^2}+\epsilon^2 \mathcal{H}_1(r,\varphi,g,L,G)+\mathcal{O}(\epsilon ^3), \end{equation}

where

(3.16)\begin{equation} {\mathcal{H}}_1(r,f,L,G)={-}\frac{1}{2}\kappa^2 r-\frac{\cos^2(f+g)}{2r}, \end{equation}

with as customary $r=a(1-e^2)/(1+e\cos f)$ and $e=\sqrt {1-(G^2/L^2)}$.

Applying theorem 2.3, we obtain the following result about the existence of symmetric periodic solutions of second kind for the anisotropic Seeliger problem (3.14).

Theorem 3.3 Fix the energy level ${\mathcal {H}}_{0}=- \frac {1}{2 L_0^2}$ and the period $T=2\pi L_0^3$ of the elliptic Kepler solution. Consider $\kappa \in \mathbb {R}^+$ and $\mu _0\in \mathbb {R}^+$ (resp. $\mu _0\in \mathbb {R}^{-}$). If $L_0> \sqrt [4]{\mu _0}/\sqrt {2\kappa }$ ( resp. $L_0> \sqrt [4]{-\mu _0}/\sqrt {2\kappa }$), then for the $2$-DOF anisotropic Seeliger's Hamiltonian (3.14) for all positive $\epsilon$ sufficiently small, there exist at least two 1-parameter or 2-parameters families of initial conditions such that each of them gives us a second-kind $S_1$ (resp. $S_2$) symmetric periodic solutions. Moreover, the periodic solutions are unstable.

Proof. We continue maintaining the notation of theorem 2.3. First, we eliminate the dependence of $\mathcal {H}_1$ given in (3.16) in the variables $r$ and $f$ using the auxiliary expressions (2.22). We infer that

(3.17)\begin{align} {\mathcal{H}}_{1}(E,g,L,G)& =\frac{1}{8 a (e \cos E-1)^3}\left[\vphantom{\left. -2 e \sqrt{1-e^2} \mu_0 \cos (\text{Ec}+2 g)\right]}2 (e \cos E-1)^2 \left(a^2 \left(e^2+2\right) \kappa ^2\right.\right.\nonumber\\ & \quad \left.\left.+\,a^2 e \kappa ^2 (e \cos (2 E)-4 \cos E)+\mu_0\right)\right.\nonumber\\ & \quad+2 e \sqrt{1-e^2} \mu_0\cos (E-2 g)-\sqrt{1-e^2} \mu_0\cos (2 (E-g))\nonumber\\ & \quad +\sqrt{1-e^2} \mu_0 \cos (2 (E+g))\nonumber\\ & \quad\left.+\,\mu_0 \cos (2 g) \left(-\left(e^2-2\right) \cos (2 E)+3 e^2-4 e \cos E\right)\right.\nonumber\\ & \quad\left. -\,2 e \sqrt{1-e^2} \mu_0 \cos (\text{Ec}+2 g)\right]. \end{align}

Next, for the calculations involving the derivative $\frac {\partial {\mathcal {H}}_{1}}{\partial G}$, we use expression (2.24) and the chain rule. We obtain

(3.18)\begin{align} \displaystyle\frac{\partial {\mathcal{H}}_{1}}{\partial G}& =\displaystyle \frac{G}{4 a^2 e (e \cos E-1)^5}\left[\frac{3}{2}\cos E(e \cos E-1)\left(2(e \cos E-1)^2\left(a^2 \left(e^2+2\right) \kappa ^2+\mu_0+a^2 e \kappa ^2\right.\right.\right.\nonumber\\ & \quad \times e \cos (2 E)-4 \cos E) +\mu_0 \cos (2 g) \left(-\left(e^2-2\right) \cos (2 E)+3 e^2-4 e \cos E\right)\nonumber\\ & \quad\left.+\, \mu_0 \sqrt{1-e^2}\left[ e \cos (E-2 g)-\cos (2 (E-g))+ \cos (2 (E+g))-2 e \cos (E+2 g)\right]\right)\nonumber\\ & \quad+\sin E\left(\frac{1}{2}\left({-}2 e (e \cos E-1)^2 \left(a^2 \left(e^2+2\right) \kappa ^2+a^2 e \kappa ^2 (e \cos (2 E)-4 \cos E)-\mu_0\right)\right.\right.\nonumber\\ & \quad+\left.\,\mu_0 \cos (2 g) \left(\left(8-12 e^2\right) \cos E+e \left(\left(2 e^2-1\right) \cos (2 E)+8 e^2-5\right)\right)\right)\sin E+3 e \left(e^2-1\right)\mu_0\nonumber\\ & \quad\times\sqrt{1-e^2} \mu_0\cos (E) \sin (2 g) \left(2 \left(e^2+1\right) \cos E-e (\cos (2 E)+3)\right)+\sin ^3(E) \cos (2 g)\nonumber\\ & \quad\left.-\,2 \sqrt{1-e^2} \mu_0 \sin ^2(E) \sin (2 g) \left({-}3 e^2+2 e \cos E+1\right)\right)-(e \cos E-1)^2\nonumber\\ & \quad\times \left(8 a^2 e^3 \kappa ^2 \cos ^4E-24 a^2 e^2 \kappa ^2 \cos ^3E+2 e \cos ^2E \left(12 a^2 \kappa ^2+\mu_0\right)\right.\nonumber\\ & \quad \left. -e\mu_0(\cos (2E)-3) \cos (2 g)\right.\nonumber\\ & \quad+\,\left.\left.2 \cos E \left({-}4 a^2 \kappa ^2+\mu_0 \left(\frac{e \sin E \sin (2 g)}{\sqrt{1-e^2}}+\sin ^2(g)-\cos ^2(g)\right)-\mu_0\right)\right.\right.\nonumber\\ & \quad\left.\left. +\frac{2 \left(1-2 e^2\right) \mu_0 \sin E \sin (2 g)}{\sqrt{1-e^2}}\right)\right]. \end{align}

Next, we evaluate (3.18) on the solution $\varphi _{0}(\tau, {\bf Y}^{j,k})$. After some simplifications, we arrive at

(3.19)\begin{align} \frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{1,k}))& = \displaystyle\frac{G_0 ({-}e + \cos E)}{16 a^2 e ( e \cos E-1)^5}\left[ 4 e \cos E \left(\mu_0-2 a^2 \left(3 e^2+4\right) \kappa ^2\right)\right.\nonumber\\ & \quad \left. +\,a^2 \left(3 \left(e^2+8\right) e^2+8\right) \kappa ^2\right.\nonumber\\ & \quad +4 (3 e^2-5) \mu_0+4 \cos (2 E) \left(e^2 \left(a^2 \left(e^2+6\right) \kappa ^2-3 \mu_0\right)+5 \mu_0\right)\nonumber\\ & \quad\left.-\,4 e \cos (3 E) \left(2 a^2 e^2 \kappa ^2+\mu_0\right)+a^2 e^4 \kappa ^2 \cos (4E)\right], \end{align}

(3.20)\begin{align} \frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{2,k}))& = \displaystyle\frac{G_0 ({-}e + \cos E)}{16 a^2 e (e \cos E-1)^5}\left[4 e \cos E \left(3\mu_0-2 a^2 \left(3 e^2+4\right) \kappa ^2\right)\right.\nonumber\\ & \quad \left. +\,a^2 \left(3 \left(e^2+8\right) e^2+8\right) \kappa ^2\right.\nonumber\\ & \quad+ 4 \left(3-4 e^2\right) \mu_0+4 \cos (2 E) \left(e^2 \left(a^2 \left(e^2+6\right) \kappa ^2+2 \mu_0\right)-5 \mu_0\right)\nonumber\\ & \quad\left.-\,4 e \cos (3 E) \left(2 a^2 e^2 \kappa ^2-\mu_0\right)+a^2 e^4 \kappa ^2 \cos (4 E)\right], \end{align}

where $a^2=L$, $G=G_0+\delta G$ and $L=L_0+\delta L$. Note that for $0< e<1$ the functions in (3.19)–(3.20) are continuous and differentiable with respect to the variable $E$ and it can be integrated on $E$. Next, for the integration of the previous two equations, we use the change (2.25). Therefore, after integration of (3.19)–(3.20) in the new variable $E$ and some simplifications, we get

(3.21)\begin{align} g^{(1)}(T/2,{\bf Y}_0^{j,k})& = \int_{0}^{\pi}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))|_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}L_{0}^{3}(1-e_0\cos E)\,{\rm d}E\nonumber\\ & =\frac{1}{2} \pi G_0 L_0^3 \left(\kappa ^2-({-}1)^{j}\frac{\left(e_0^2+2 \sqrt{1-e_0^2}-2\right) \mu_0}{a^2 e_0^4 \sqrt{1-e_0^2}}\right). \end{align}

Replacing $e_0=\sqrt {L_0^2-G_0^2}/L_0$ on equation (3.21), we obtain

(3.22)\begin{align} & g^{(1)}(T/2,{\bf Y}_0^{j,k})\nonumber\\ & \quad =\frac{\pi G_0 L_0^2 \left(G_0^4 |G_0|\kappa ^2 L_0 -2 G_0^3 \kappa ^2 L_0^3+|G_0| \kappa ^2 L_0^5-({-}1)^j\mu_0( G_0^2+L_0^2-2L_0|G_0|)\right)}{2 |G_0| (G_0-L_0)^2 (G_0+L_0)^2}. \end{align}

Now, we must search the solutions $G_0$ of the equation $g^{(1)}(T/2,{\bf Y}_0^{j,k})=0$, separately, in the two domains

(3.23)\begin{equation} D_1=\{(L_0,G_0) \in \mathbb{R}^2; 0< G_0< L_0\}, \quad D_2=\{(L_0,G_0) \in \mathbb{R}^2\in \mathbb{R}; -L_0< G_0<0\}. \end{equation}

First, in the domain $D_1$( $G_0>0$) we need to solve the equation

(3.24)\begin{equation} g^{(1)}(T/2,{\bf Y}_0^{j,k})= \frac{\pi L_0^2 \left[L_0 \kappa^2 G_0^3+2L_0^2\kappa^2G_0^2+L_0^3\kappa ^2G_0+({-}1)^j\mu_0\right]}{2 (G_0+L_0)^2}=0. \end{equation}

The solutions of equation (3.24) are the roots of the polynomial equation

(3.25)\begin{equation} p_j(G_0)= L_0 \kappa^2 G_0^3+2L_0^2\kappa^2G^2+L_0^3\kappa ^2G_0+({-}1)^j\mu_0, \end{equation}

in the interval $I_p=[0,L_0]$. Initially, we consider the case $j=1$. If $\mu _0<0$, then equation (3.25) has no solution. A simple inspection shows that

\[ p_1'(G_0)=L_0\kappa^2(3 G_0^2+4 L_0 G_0+ L_0^2) >0, \text{ for all } G_0\in I_p. \]

Moreover, if $\mu _0>0$ and $L_0 > \mu _0^{1/4}/\sqrt {2\kappa }$, then $p_1(0)=-\mu _0<0$ and $p_1(L_0)=-\mu _0 + 4 \kappa ^2 L_0^4>0$. Thus, for $j=1$, $\mu _0>0$ and $L_0 > \mu _0^{1/4}/\sqrt {2\kappa }$, there exists only one solution $G_{0,1}=G_{0,1}(\kappa,L_0,\mu _0)\in (0,L_0)$ of the polynomial equation $p_1(G_0)$.

On the other, for the case $j=2$, $\mu _0<0$ and $L_0 > (-\mu _0)^{1/4}/\sqrt {2\kappa }$, there exists only one solution $G_{0,2}=G_{0,2}(\kappa,L_0,\mu _0)\in (0,L_0)$ of the polynomial equation $p_2(G_0)$.

Next, we verify the condition of non-degeneracy (b) of theorem 2.3. Again, we use the auxiliary expression (2.22) and the chain rule for the calculus involving the derivative of equations (3.19)–(3.20) with respect to $\delta G$. After that, we evaluate $\frac {\partial }{\partial \delta G}\left (\frac {\partial \mathcal {H}_1}{\partial G}\right )$ on the solution $\varphi _{0}(t, {\bf Y}_{0}^{j,k}))$. Making the change of variable $t=L_{0}^{3}(E-e_0\sin E)$, after the integration $\frac {\partial ^2 \mathcal {H}_1}{\partial G\partial \delta G}(\varphi _{0}(t, {\bf Y}_0^{j,k}))$ and some simplification, we arrive to

(3.26)\begin{align} & \frac{\partial g^{(1)}}{\partial \delta G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}\nonumber\\ & \quad =\int_{0}^{\pi}\frac{\partial }{\partial \delta G}\left(\frac{\partial \mathcal{H}_1}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\right)\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}L_0^3(1-e_0\cos E)\,{\rm d}E\nonumber\\ & =\pi \frac{L_0^3}{2 a^3 e_0^6} \left(a^3 e_0^6 \kappa ^2+({-}1)^{j+1}\frac{a \left(e_0^2+2 \sqrt{1-e_0^2}-2\right) e_0^2 \mu_0}{\sqrt{1-e_0^2}}({-}1)^{j+1}\right.\nonumber\\ & \quad\left. \frac{ \left(3 e_0^4+4 \left(2 \sqrt{1-e_0^2}-3\right) e_0^2-8 \sqrt{1-e_0^2}+8\right) G_0^2\mu_0}{\left(1-e_0^2\right)^{3/2}}\right) \end{align}

Replacing $a=L_0^2$ and $e_0=\sqrt {L_0^2-G_0^2}/L_0$ on equation (3.26), we obtain

(3.27)\begin{align} & \frac{\partial g^{(1)}}{\partial \delta G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}\nonumber\\ & \quad =\int_{0}^{\pi}\frac{\partial }{\partial \delta G}\left(\frac{\partial \mathcal{H}_1}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\right)\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}L_0^3(1-e_0\cos E)\,{\rm d}E\nonumber\\ & =\frac{\pi L_0^3}{2 \left(L_0^2-G_0^2\right)^3} \left[({-}1)^{j+1}\frac{\mu_0 \left(G_0^2-L_0^2\right) \left(G_0^2+L_0^2-2 L_0 | G_0|\right)}{|G_0| L_0}\right.\nonumber\\ & \left.\quad +\,({-}1)^j\frac{ \mu_0 \left(6 G_0^2 L_0^2-8 G_0^2 L_0 | G_0| +3 G_0^4-L_0^4\right)}{|G_0| L_0}\right.\nonumber\\ & \quad\left.-\,\kappa ^2 \left(G_0^2-L_0^2\right)^3\right] \end{align}

Considering $G_0>0$ in (3.27) we infer that

(3.28)\begin{equation} \frac{\partial g^{(1)}}{\partial \delta G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}=\displaystyle\frac{\pi L_0^2 \left[({-}1)^{j-1}2\mu_0+\kappa^2 L_0 (G_0+L_0)^3\right]}{2(G_0+L_0)^3}. \end{equation}

Now, if $\mu _0\!>\!0$ (resp. $\mu _0\!<\!0)$, then $\frac {\partial g^{(1)}}{\partial \delta G}(\varphi _{0}(\tau, {\bf Y}_0^{1,k}))\!>\!0$ (resp. $\frac {\partial g^{(1)}}{\partial \delta G}(\varphi _{0}(\tau, {\bf Y}_0^{2,k}))\!>\!0$) for all $G_0\in (0,L_0)$. Therefore, by theorem 2.3, we obtain the existence of a $1$-parameter (on $\epsilon$) family (and a $2$-parameter on $\epsilon$ and $\delta L$) of initial conditions [as in (2.15) and (2.16)] such that each of them gives rise to $S_j$-symmetric (prograde) periodic solution.

Second, for the study of the retrograde solutions ($G_0<0$), for each $j=1,2$, we need to search solutions of (3.22) in the domain $D_2$. Equation (3.22) for $G_0<0$ assumes the form

(3.29)\begin{equation} g^{(1)}(T/2,{\bf Y}_0^{j,k})= \frac{\pi L_0^2 \left[L_0 \kappa^2 G_0^3-2L_0^2\kappa^2G_0^2+L_0^3\kappa ^2G_0+({-}1)^{j-1}\mu_0\right]}{2 (G_0-L_0)^2}. \end{equation}

The negative solutions on the variable $G_0$ of $g^{(1)}(T/2,{\bf Y}_0^{j,k})=0$ are the solutions of the polynomial equation

(3.30)\begin{equation} q_j(G_0)=L_0 \kappa^2 G_0^3-2L_0^2\kappa^2G^2+L_0^3\kappa ^2G_0+({-}1)^{j-1}\mu_0=0. \end{equation}

It is easy to check that $q_j(G_0)=-p_j(-G_0)$, $-L_0< G_0<0$ and $j=1,2$. Thus, for each $j=1,2$, there exists a unique negative solution

\[ \overline{G}_{0,j}=\overline{G}_{0,j}(\kappa,L_0,\mu_0)={-}G_{0,j}(\kappa,L_0,\mu_0), \]

defined in the interval $(-L_0,0)$, since for $j=1$, $L_0 > \mu _0^{1/4}/\sqrt {2\kappa }$ and for $j=2$, $L_0 > (-\mu _0)^{1/4}/\sqrt {2\kappa }$.

Now, we analyse the condition (b) given in theorem 2.3. To calculate $\frac {\partial }{\partial \delta G}\left (\frac {\partial \mathcal {H}_1}{\partial G}\right )(\varphi _{0}(t, {\bf Y}_{0}^{j,k}))$, we proceed as in the previous case and after some manipulation we get

(3.31)\begin{align} \frac{\partial g^{(1)}}{\partial \delta G}\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}& = \displaystyle\int_{0}^{\pi}\frac{\partial }{\partial \delta G}\left(\frac{\partial \mathcal{H}_1}{\partial G}\right)(\varphi_{0}(t, {\bf Y}^{j,k}))\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}L_0^3(1-e_0\cos E)\,{\rm d}E\nonumber\\ & = \frac{\pi L_0^2 \left[({-}1)^{j-1}2\mu_0+\kappa^2 L_0 (L_0-{G_0})^3\right]}{2(L_0-G_0)^3}. \end{align}

If $\mu _0>0$ (resp. $\mu _0<0)$, then $\frac {\partial g^{(1)}}{\partial \delta G}(\varphi _{0}(t, {\bf Y}_0^{1,k}))>0$ (resp. $\frac {\partial g^{(1)}}{\partial \delta G}(\varphi _{0}(t, {\bf Y}_0^{2,k}))>0$) for all $G_0\in (-L_0,0)$. Thus, the conditions (a) and (b) of theorem 2.3 are satisfied and it is guaranteed of the existence of a $1$-parameter (on $\epsilon$) (resp. a $2$-parameter on $\epsilon$ and $\delta L$) family of initial conditions such that each of them gives us $S_j$-symmetric retrograde periodic solution of the Hamiltonian (3.14).

Now, we analyse the type of stability of the previous symmetric periodic solutions. After some algebraic manipulation, the averaging function $\overline {H}$ given in (2.10) for the Hamiltonian problem (3.14), in Delaunay variables, has the form

(3.32)\begin{align} \overline{H}(X_1,X_2)& =\frac{-\pi}{2} L_0 \left[\mu_0\frac{\cos (2 (g_0^{j,k}+X_1)) \left(( G_0+X_2)^2-2L_0 | G_0+X_2| +L_0^2\right)}{L_0^2-( G_0+X_2)^2}\right.\nonumber\\ & \left.\quad+\,\kappa ^2L_0^2 \left(3 L_0^2-( G_0+X_2)^2\right)+\mu_0\vphantom{\left[\mu_0\frac{\cos (2 (g_0^{j,k}+X_1)) \left(( G_0+X_2)^2-2L_0 | G_0+X_2| +L_0^2\right)}{L_0^2-( G_0+X_2)^2}\right.}\right]. \end{align}

First, we consider equation (3.32) for prograde solutions ($G_0>0$). Therefore, the matrix $A$ defined in (2.11) takes the form

(3.33)\begin{align} A_P^{(j)} & = \mathbb{J}\left( \frac{\partial^2 \bar{H}}{\partial X_k \partial X_l} \right)_{X=0}\nonumber\\ & \quad =\left( \begin{array}{cc} 0 & \displaystyle L_0^2\pi\left(L_0\kappa^2+\dfrac{2({-}1)^{j-1}\mu_0}{(G_0+L_0)^3}\right) \\ \dfrac{2 \pi L_0 ({-}1)^{j-1}\mu_0 (L_0-G_0)}{G_0+ L_0} & 0 \\ \end{array}\right), \end{align}

for $j=1,2$, whose eigenvalues are

\[ \lambda_{1,2}^{(j)}={\pm} \frac{ \sqrt{2} \pi L_0^{3/2} \sqrt{({-}1)^{j-1}\mu_0(L_0-G_0)} \sqrt{\kappa^2 L_0 (G_0+L_0)^3+2({-}1)^{j-1}\mu_0}}{(G_0+L_0)^2}. \]

Since $0< G_0< L_0$, then $\kappa ^2 L_0 (G_0+L_0)^3+2(-1)^{j-1}\mu _0>0$ and $(-1)^{j-1}\mu _0 (L_0-G_0) >0$ for $j=1$ and $\mu _0>0$ or $j=2$ and $\mu _0<0$. So, the prograde periodic solutions have multipliers characteristic

\[ 1,1,1+\epsilon \lambda_1^{(j)}+ \mathcal{O}(\epsilon^2),1+\epsilon \lambda_2^{(j)}+ \mathcal{O}(\epsilon^2), \]

and therefore these $S_j$-symmetric periodic solutions are unstable.

On the other, considering equation (3.32) with $G_0<0$ implies that the matrix $A$ takes the form

(3.34)\begin{align} A_R^{(j)}& = \mathbb{J}\left( \frac{\partial^2 \bar{H}}{\partial X_k \partial X_l}\right)_{X=0}\nonumber\\ & \quad =\left( \begin{array}{cc} 0 & \displaystyle L_0^2\pi\left(L_0\kappa^2+\dfrac{2({-}1)^{j-1}\mu_0}{(L_0-G_0)^3}\right) \\ \dfrac{2 \pi L_0 ({-}1)^{j-1}\mu_0 (L_0+G_0)}{L_0-G_0} & 0 \end{array}\right), \end{align}

for $j=1,2$. The eigenvalues of the matrix $A_R^{(j)}$ are given by

\[ \rho_{1,2}^{(j)}={\pm} \frac{ \sqrt{2} \pi L_0^{3/2} \sqrt{({-}1)^{j-1}\mu_0(L_0+G_0)} \sqrt{\kappa ^2 L_0 (L_0-G_0)^3+2({-}1)^{j-1}\mu_0}}{(G_0-L_0)^2}. \]

Since $-L_0< G_0<0$, then $\kappa ^2 L_0 (L_0-G_0)^3+2(-1)^{j-1}\mu _0>0$ and $(-1)^{j-1} \mu _0(L_0+G_0) >0$ for $j=1$ and $\mu _0>0$ or $j=2$ and $\mu _0<0$. So, the eigenvalues $\rho _{1,2}^{(j)}$ are reals and the retrograde periodic solutions have characteristic multipliers

\[ 1,1,1+ \epsilon \rho_{1,2}^{(j)}+\mathcal{O}(\epsilon^2),1- \epsilon \rho_{1,2}^{(j)}+\mathcal{O}(\epsilon^2). \]

Thus, we conclude that the retrograde symmetric periodic solutions are unstable.

3.3 The planar generalized Størmer problem

This problem consists in the study of the dynamics of a charged particle around rotating magnetic planets. More specifically, the generalized Størmer problem describes the dynamics of a dust particle of mass $m$ and charge $q$ orbiting a rotating magnetic planet of mass $M$. The magnetic field of the planet is supposed to be a perfect magnetic dipole of strength aligned along the north–south poles of the planet (see [17–19] for more details on the formulation of this problem). Moreover, the planet magnetosphere is taken as a rigid conducting plasma which rotates with the same angular velocity $\Omega$ as the planet, in such a way that the charge $q$ is subject to a co-rotational electric field. Furthermore, the gravitational interaction in this model takes into account the non-sphericity of the planet which is given by means of the so-called $J_2$ term. See [19] for a complete discussion of the Størmer problem with $J_2$ effect. We will consider the planar generalized Størmer problem which is given through the following two-degree-of-freedom Hamiltonian

(3.35)\begin{align} H(x,y,p_x,p_y)& =\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{1}{\sqrt{x^2+y^2}}-\frac{\delta (xp_y-yp_x)}{(x^2+y^2)^{3/2}}+ \frac{\delta\beta}{\sqrt{x^2+y^2}}\nonumber\\ & \quad-\frac{J_2}{2(x^2+y^2)^{3/2}} +\frac{\delta^2}{2(x^2+y^2)^2}. \end{align}

It is easy to check that the Hamiltonian function (3.35) is invariant under the two symmetries $S_1$ and $S_2$. The planar generalized Størmer problem (3.35) depends also on three external parameters, namely, $\delta, \beta$ and $J_2$. The parameter $\delta$ indicates the ratio between the magnetic and the Keplerian interaction (i.e., the charge–mass ratio $q/m$ of the particle). The parameter $\beta$ is the ratio between the electrostatic and the Keplerian interactions (i.e., the ratio $\Omega /w_k$, where $w_k=\sqrt {M/\mathcal {R}}$ and $\mathcal {R}$ is the equatorial radius of the planet). Finally, $J_2$ is the oblateness of the planet taken into consideration.

For our purpose, we introduce the small parameter $\epsilon$ by means of the following relations $\delta = \epsilon b$ and $J_2=\alpha \epsilon$ with $b, \alpha$ real numbers. So, the Hamiltonian function (3.35) becomes

(3.36)\begin{align} H& =\frac{1}{2}(p_{x}^2+p_{y}^2)-\frac{1}{\sqrt{x^2+y^2}}\nonumber\\& \quad +\epsilon\left[\frac{-b (xp_y-yp_x)}{(x^2+y^2)^{3/2}}+ \frac{b\beta}{\sqrt{x^2+y^2}}-\frac{\alpha}{2(x^2+y^2)^{3/2}}\right]+\mathcal{O}(\epsilon^2).\end{align}

The next step is to express (3.36) in the mixed coordinates involving polar and Delaunay elements. Then, we arrive at

(3.37)\begin{equation} \mathcal{H}={-}\frac{1}{2L^2}+\mathcal{H}_1(r,\varphi,g,L,G)+\mathcal{O}(\epsilon ^2), \end{equation}

where $(r,\theta )$ are the classical polar coordinates, $\theta =f+g$ and $f$ is the true anomaly. The perturbed function $\mathcal {H}_1$ is given by

(3.38)\begin{equation} {\mathcal{H}}_1={-}\frac{b G}{r^3}+\frac{b\beta}{r}-\frac{\alpha}{2 r^3}={-}\frac{-2 b \beta L^2 (1- e \cos E)^2+\alpha +2 b G}{2 L^2 (1- e\cos E)^3}, \end{equation}

where $r=a(1-e^2)/(1+e\cos f)$ and $e=\sqrt {1-(G^2/L^2)}$ is the eccentricity of the unperturbed elliptic orbit.

To obtain symmetric periodic solutions, we must verify the conditions (a) and (b) in theorem 2.3. We have the following result for the existence of the symmetric periodic solutions of second kind for the generalized Størmer problem (3.36).

Theorem 3.4 Fix the energy level ${\mathcal {H}}_{0}= - \frac {1}{2 L_0^2}$ and the period $T=2\pi L_0^3$ of the elliptic Kepler solution. Given $\alpha$ and $b$ non null real constants, for the $2$-DOF generalized Størmer problem (3.36) for all $\epsilon$ and $\delta L$ positive and sufficiently small, the following statement holds:

If $\alpha b<0$ (resp. $\alpha b> 0$) and $L_0> \frac {3}{4}\left |-\frac {\alpha }{b}\right |$, then there exist two $1$-parameter ($\epsilon$) families and two $2$-parameter ($\epsilon$ and $\delta L$) families of initial conditions such that each of them gives rise to a prograde (resp. retrograde) second-kind $S_i$-symmetric periodic solution.

The $S_1$-symmetric periodic solutions are obtained as continuation of the elliptic Keplerian solution with initial conditions

\[ Y_{0}^{1,k}=\left(0,(2k+1)\frac{\pi}{2},L_0, G_0\right), \quad k=0,1, \]

where $G_0=-\frac {3\alpha }{4b}$. On the other, the $S_2$-symmetric periodic solutions are obtained as continuation of the elliptic Keplerian solution with initial condition

\[ Y_{0}^{2,k}=\left(0,k\pi,L_0, G_0\right),\quad k=0,1. \]

Moreover, the families of symmetric periodic solutions generated by the $1$-parameter families of initial conditions have fixed period $T=2\pi L_{0}^{3}$ and those generated by the $2$-parameter families of initial conditions have period $\overline {T}=T(1- \epsilon T^*)+\mathcal {O}(\epsilon ^2)$ where $T^*=\frac {4}{9} \pi b L_0 (9 \beta -\frac {8 b^2 L_0}{\alpha ^2})$. All these symmetric periodic solutions are close to elliptic Keplerian solutions with eccentricity $e_0=\sqrt {16b^2L_0^2-9\alpha ^2}/4bL_0$.

Proof. Maintaining the notation of § 2, we will verify the hypotheses (a) and (b) of theorem 2.3. For the condition (a), we need to find the solutions of

(3.39)\begin{equation} g^{(1)}(T/2,{\bf Y}^{j,k})\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}=\int_{0}^{T/2}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}^{j,k}))\Big|_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}d\tau=0, \end{equation}

for $j=1$ and $j=2$. Observe that the perturbed function $\mathcal {H}_1$ in (3.38) does not depend on the variable $g$, and therefore, we have $\frac {\partial {\mathcal {H}}_{1}}{\partial G}(\varphi _{0}(\tau, {\bf Y}_{0}^{1,k}))= \frac {\partial {\mathcal {H}}_{1}}{\partial G}(\varphi _{0}(\tau, {\bf Y}_{0}^{2,k}))$ for $k=1,2$. Using expression (2.24) for calculations involving the derivative $\frac {\partial {\mathcal {H}}_{1}}{\partial G}$, we arrive at

(3.40)\begin{align} \displaystyle \frac{\partial {\mathcal{H}}_{1}}{\partial G}& = \displaystyle \frac{1}{4 a^4 e (e \cos E-1)^5} \left[\cos E \left(a^2 b \beta \left(11 e^2+4\right) G-8 a b e^2-6 G (\alpha +2 b G)\right)\right.\nonumber\\ & \quad\left. +\,e\left(4ab+2 a b e^2 -2 a^2 b \beta e^2 G-8 a^2 b \beta G-2 a b \cos (2 E) \left(e^2 (a \beta G-1)\right.\right.\right.\nonumber\\ & \quad\left.\left.\left. +\,2 a \beta G\right)+12 b G^2+6 \alpha G+a^2 b \beta e G \cos (3 E)\right)\right]. \end{align}

where $a=L^2$. Next, we evaluate (3.40) on the solution $\varphi _{0}(t, {\bf Y}_{0}^{j,k})$. It follows that

(3.41)\begin{align} g^{(1)}(T/2,{\bf Y}_0^{j,k})& =\displaystyle\int_{0}^{T/2}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}_0^{j,k}))\,{\rm d}\tau\nonumber\\ & = \int_{0}^{\pi}\frac{\partial {\mathcal{H}}_{1}}{\partial G}(\varphi_{0}(\tau, {\bf Y}_0^{j,k}))L_0^3(1-e_0\cos E)\,{\rm d}E, \end{align}

Therefore, after integration of (3.41) we arrive at

(3.42)\begin{equation} g^{(1)}(T/2,{\bf Y}_0^{j,k})=\frac{ L_0^3 \left(2 a b \left( e_0^2-1\right)+3 G_0 (\alpha +2 b G_0)\right)}{2 a^4 \left(1-e_0^2\right)^{5/2}}=\frac{\pi (3 \alpha +4 b G_0)}{2 G_{0}^3 | G_0| }. \end{equation}

We must find the solutions $G_0$ of the $g^{(1)}(T/2,{\bf Y}_0^{j,k})=0$, in the two domains given in (3.23). Observe that in both domains, the solution of the equation $g^{(1)}(T/2,{\bf Y}_0^{j,k})=0$ is given by

(3.43)\begin{equation} G_0={-}\frac{3\alpha}{4b}. \end{equation}

Thus, in the region $D_1$ (prograde solutions), we need to impose that $\alpha b<0$ and $L_0>-\frac {3 \alpha }{4 b}$. While, in the region $D_2$ (retrograde solutions), we must have $\alpha b>0$ and $L_0>\frac {3 \alpha }{4 b}$.

Now, we verify the condition (b) of theorem 2.3. To compute the partial derivative $\frac {\partial }{ \partial \delta G}\left (\frac {\partial \mathcal {H}_1}{ \partial G}\right )$, we use the chain rule and expressions (2.24) to differentiate (3.40) with respect to $\delta G$. Thus, again using the change of variables $t=L_{0}^{3}(E-e_0\sin E)$, after integration, we obtain

(3.44)\begin{equation} \frac{\partial g^{(1)}}{\partial \delta G}\Big |_{{\bf Y}^{j,k}\!=\!{\bf Y}_{0}^{j,k}}=\int_{0}^{\pi}\frac{\partial }{ \partial \delta G}\left(\frac{\partial \mathcal{H}_1}{ \partial G}\right)\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}L_0^3(1\!-\!e_0\cos E)\,{\rm d}E\!={-}\frac{6 \pi (\alpha +b G_0)}{G_0^4 |G_0|}. \end{equation}

Evaluating expression (3.44) on the solutions $G_0=-\frac {3 \alpha }{4 b}$, we get

\[ \frac{\partial g^{(1)}}{\partial \delta G}\Big |_{{\bf Y}^{j,k}={\bf Y}_{0}^{j,k}}=\frac{512 \pi b^5}{81 \alpha ^4}\neq 0, \]

since $b\neq 0$. Therefore, we have verified the conditions given in (2.12). Thus, we conclude that every Keplerian elliptic orbit with initial condition

\[ {\bf Y}_{0}^{j,k}= \left(0, g_{0}^{j,k}, L_0,\pm G_0\right),\quad G_0={-}\frac{3\alpha}{4b}, \quad k=1,2, \]

can be continued to a $S_j$-symmetric periodic solution of the generalized Størmer problem. Therefore, we conclude the proof.

We point out that in the previous theorem, we cannot give information about the linear stability of the symmetric periodic solutions, because the matrix $A$ (2.11) is given by $\left (\begin {array}{cc} 0 & * \\ 0 & 0 \\ \end {array}\right )$, so their eigenvalues are all null.

Remark 3.5 A simple calculation shows that the averaged function of (3.38) in Delaunay variables is

\[ \overline{\mathcal{H}}=\frac{-\alpha +2 b \beta L | G| ^3-2 b G}{2 L^3 | G| ^3}. \]

Moreover, the averaged system associated with the Hamiltonian (3.36), for prograde solutions is given by

(3.45)\begin{equation} \begin{aligned} \frac{{\rm d} G}{{\rm d} \ell} & ={-}\epsilon \frac{1}{({-}2h^*)^{3/2}}\frac{3 \alpha +4 b \text{G}}{2 \text{G}^4 \text{L}^3},\\ \frac{{\rm d} g}{{\rm d} \ell} & =0. \end{aligned} \end{equation}

Therefore, system (3.45) has no non-degenerate critical points. The same conclusion is valid for retrograde solutions, and Reeb's Theorem or Averaging theory does not give us information about the existence of periodic solutions.

4. Concluding remarks

In this work, we have considered the problem of existence of periodic solutions of symmetric Hamiltonian systems which are perturbation of the integrable Kepler problem with 2-DOF using an analytic approach. The method of analysis depends on the appropriated use of Delaunay coordinates, because they permit us to get the first approximation of the solutions of the full Hamiltonian system. This is achieved as solution of a variational system. Next, we give the sufficient conditions for the existence of second-kind symmetric periodic solutions for problem (1.1) as continuation of an elliptic Kepler solution (the so-called continuation Poincaré's method). Moreover, we get an estimate of the characteristic multipliers of the symmetric periodic solutions, thus allowing us to determine the type of stability of such solutions.

We complement our study considering the connection between the solutions obtained in this paper and the procedure that allows to get symmetric periodic solutions applying the averaging theory for Hamiltonian systems and symplectic reduction. For this last approach, the main result can be found in [36], here the authors derive a method to get symmetric periodic solutions from a more general class of periodic solutions that are obtained from the analysis of relative equilibria after performing the process of averaging and reduction. In this work, we determine when an elliptic Keplerian solution that can be continued by theorem 2.3 to a symmetric periodic solution of the perturbed problem (1.1) can also be continued from the process of the averaging (Reeb's Theorem). As we saw in § 3.3, some degenerate elliptic Keplerian solutions (associated with a degenerate critical point in the sense of Reeb's theorem) can be continued by theorem 2.3. In addition, if a Keplerian elliptic solution can be continued by theorem 2.3 and by Reeb's Theorem to a periodic solution of the full problem, then we showed that these periodic solutions have the same linear stability.

As applications of our theoretical results, we study the existence of periodic solutions of three different problems: the perturbed hydrogen atom with stark and quadratic Zeeman effect, for the anisotropic Seeligers two-body problem and to the planar generalized Størmer problem. We proceeded with providing care when carrying out our analysis, providing many but necessary details to apply theorem 2.3. After checking the literature on the subject, we emphasize that our result stated on theorems 3.1, 3.3 and 3.4 are new.

In a future work, we intend to study the existence of second-kind symmetric periodic solutions for spatial perturbed Kepler problems. In addition, we want to establish the connection between the solutions obtained by the classical analytical continuation method of Poincaré and the procedure to obtain symmetric periodic solutions applying averaging theory of Hamiltonian systems and symplectic reduction for three degrees of freedom.