Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 63, pp. 1–42.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

PERIODIC ORBITS OF THE SPATIAL ANISOTROPIC KEPLER

PROBLEM WITH ANISOTROPIC PERTURBATIONS

MENGYUAN LI, QIHUAI LIU

Abstract. In this article, we study the periodic orbits of the spatial anisotropic
Kepler problem with anisotropic perturbations on each negative energy sur-

face, where the perturbations are homogeneous functions of arbitrary integer

degree p. By choosing the different ranges of a parameter β, we show that
there exist at least 6 periodic solutions for p > 1, while there exist at least 2

periodic solutions for p ≤ 1 on each negative energy surface. The proofs of

main results are based on symplectic Delaunay coordinates, residue theorem,
and averaging theory.

1. Introduction

Celestial mechanics has stimulated the development of many branches of Math-
ematics [12, 14]. The two-body problem or Kepler problem is the basic model of
celestial mechanics, which attracts the vivid interest of many mathematicians who
have studied its classical form with different types of perturbations. The anisotropic
problem arises in many areas, and it studies the motion of particles, where the in-
teraction law is different in each direction of the space. The study of the dynamics
of the planar anisotropic celestial problem mainly includes the anisotropic Manev
problem and the anisotropic Kepler problem with perturbations, see [24, 31].

The plane anisotropic Kepler problem with perturbations is described by the
Hamiltonian

H =
1

2

(
X2 + Y 2

)
− 1√

µx2 + y2
− εβ

(µx2 + y2)p/2
, (1.1)

where β is a constant, p is a integer and ε is a small parameter.
The Hamiltonian (1.1) for β = 0 and µ = 1 corresponds to the classical Kepler

problem. In this respect, quite mature work has been done. See for example the
books [5, 38] for a detailed introduction. If β 6= 0 and µ = 1, we have the Kepler
problem with perturbations which are symmetric with respect to the origin. Vidal
[46] proved that each circular solution of the unperturbed problem gives rise to a
periodic solution of the perturbed system.

2010 Mathematics Subject Classification. 37N05, 37G15, 70F05.
Key words and phrases. Periodic orbit; averaging theory; residue theorem;

spatial anisotropic Kepler problem.
c©2021. This work is licensed under a CC BY 4.0 license.

Submitted April 12, 2020. Published July 8, 2021.

1



2 M. LI, Q. LIU EJDE-2021/63

For β = 0 and µ 6= 1, the Hamiltonian (1.1) becomes

H =
1

2

(
X2 + Y 2

)
− 1√

µx2 + y2
, (1.2)

and we have the anisotropic Kepler problem which comes originally from quantum
mechanics. It has been studied by Gutzwiller to establish a link between quantum
mechanics and classical mechanics, see [21, 20, 22, 24], and these results were later
extensively extended by Devaney [8, 9, 10], Casasayas and Lliber [3]. Devaney
[8] studied the collision manifold of the anisotropic Kepler problem and obtained
that, for all µ > 1, the primary bi-collision orbits along the x-axis are transversal
heteroclinic solutions of the system. Moreover, based on the invariant boundary
(McGehee’s method) to “blown up” the singularity, he obtained a qualitative pic-
ture of the behavior of a mechanical system near a singularity, see [9, 10]. Casasayas
and Llibre exhibited phenomena such as non-integrability and chaotic behaviour,
and also surveyed the recent techniques and results from the anisotropic Kepler
problem, see [3]. By using the Lyapunov characteristic number (LCN), Contopou-
los and Harsoula studied the degree of order and chaos for the anisotropic Kepler
problem, calculate the asymptotic curves from unstable periodic orbits, and conjec-
ture that there are stable orbits of high multiplicities for larger values of µ. These
curves go through the collision orbits and continue beyond these points, see [4].

The existence of periodic orbits of the anisotropic Kepler problem is of special
interest for many mathematicians. Firstly, Gutzwiller studied numerical periodic
solutions of the plane anisotropic Kepler problem, and show that periodic orbits in
the plane anisotropic Kepler problem with mass-ratio five are found numerically up
to length five; see [23] for details. By using the McGehee coordinates and averaging
theory, Abouelmagd, Llibre and Guirao [1] proved that at every energy level the
anisotropic Kepler problem with small anisotropy has at least two periodic orbits.
These result has been generalized by Lliber and Valls in [31].

When β 6= 0, µ = 1 and p = 2, we have

H =
1

2

(
X2 + Y 2

)
− 1√

x2 + y2
− εβ

x2 + y2
(1.3)

which corresponds to the Manev problem. One of the advantages of the Manev
problem over the Keplerian is that it explains the perihelion advance of the inner
planets with the same accuracy as relativity [44]. Firstly, Delgado et al. founded an-
alytic solutions of Manev systems corresponding to (1.3) and described completely
the global flows by using the McGehee coordinates and topological methods, see [7].
They showed that if the energy constant is negative, then the orbits are generically
precessional ellipses; for zero energy, the orbits are precessional parabolas, and for
positive energy, they are precessional hyperbolas. Lliber, Teruel and Fuente [30]
characterized the global flows of Manev systems, and gave the phase portraits.
Szenkovits, Mioc and Stoica [42] tackled the Manev problem from the standpoint
of topology and pointed out the first integrals of energy and angular momentum.
We can refer to [17, 34, 35, 36, 37, 39, 40, 43] for related developments.

When β 6= 0, µ 6= 1 and p = 2, the corresponding problem is called the
anisotropic Manev problem (AMP) which is inspired by the anisotropic Kepler
problem. One of the main purpose for studying the anisotropic Manev problem, is
to further analyze the similarities between classical mechanics and quantum the-
ory. This was introduced by Craig, Diacu, Lacomba and Pérez in [6]. With the



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 3

McGehee coordinates, they found a positive-measure set of collision orbits and tack-
led capture and escape solutions in the zero-energy case. Santoprete [41] studied
the existence of periodic solutions for weak anisotropy and found that the sym-
metric periodic orbits of the Manev system are perturbed to periodic orbits in
the anisotropic problem. Moreover, Diacu and Santoprete proved that for weak
anisotropy, chaos shows up on the zero-energy manifold by using an extension of
the Poincaré-Melnikov method, see [15]. Recently, Based on the method of averag-
ing, Llibre and Yuan showed that the Kepler problem has a unique elliptic periodic
solution can be continued to the plane anisotropic Manev problem under a certain
condition, see [32] for details. We refer to the related papers [11, 16].

If β 6= 0, µ 6= 1, two kinds of models are considered. On the one hand, the plane
Kepler problem with anisotropic perturbations is given by the Hamiltonian

H =
1

2

(
X2 + Y 2

)
− 1√

x2 + y2
− εβ

(µx2 + y2)
p/2

, (1.4)

has been first studied by Diacu, Pérez-Chavela and Santoprete in [13], where µ is
near 1 and β is a constant, p is a integer and ε is a small parameter. They obtained
that for p > 2, the sets of initial conditions leading to collisions/ejections and the
one leading to escapes/captures have positive measure; for p > 2 and p 6= 3, the
flows on the zero-energy manifold are chaotic and for p = 2, the infinity manifold of
the zero-energy level has heteroclinic connections with the collision manifold. The
results are extended by Escalona-Buendia and Pérez-Chavela in [18]. On the other
hand, the plane anisotropic Kepler problem with anisotropic perturbation has been
considered by Miguel, Raquel and Juan in [33]. Based on averaging theory, by
using symplectic Delaunay coordinates they analyzed the sufficient conditions for
existence and kind of stability of periodic orbits.

The planar models in celestial mechanics are naturally extended to the spatial
cases and there also are a few results about the dynamics of the spatial anisotropic
problem. See for example[19, 29]. The spatial anisotropic Kepler problem with
anisotropic perturbations is given by the Hamiltonian

H =
1

2

(
X2 + Y 2 + Z2

)
− 1√

µ
(
x2 + y2

)
+ z2

− εβ(
µ
(
x2 + y2

)
+ z2

)p/2 , (1.5)

where β is a constant, p is a integer and ε is a small parameter.
The spatial Kepler problem and the spatial Manev problem are corresponding to

the Hamiltonian (1.5) with the parameters β = 0, µ = 1 and β 6= 0, µ = 1, p = 2,
respectively. One of the advantages of the Manev problem over the Keplerian is
that it explains the perihelion advance of the inner planets with the same accuracy
as relativity [45].

If β = 0 and µ 6= 1, we have spatial anisotropic Kepler problem. The existence of
periodic orbits is firstly studied by Guriao, Llibre and Vera in [19], where periodic
orbits of the perturbed spatial Keplerian Hamiltonians with axial symmetry, such as
Matese-Whitman Hamiltonian and spatial generalized van der Waals Hamiltonian
were considered.

If β 6= 0, µ 6= 1 and p = 2, we have spatial anisotropic Manev problem (SAMP).
Llibre and Makhlouf have proved the existence of periodic orbits on the every level
H = h < 0 by using the method of averaging in [29]. Later, the results were



4 M. LI, Q. LIU EJDE-2021/63

generalized by Alberti and Vidal in [2] by using symmetry conditions. They found
the existence of new families of periodic orbits.

In this article, inspired by the study of plane anisotropy Kepler problem with
anisotropy perturbations, we study the existence of periodic orbits on negative
energy level H = h (h < 0) for the spatial anisotropic problem with the Hamiltonian

H =
1

2

(
X2 + Y 2 + Z2

)
− 1√

µ
(
x2 + y2

)
+ z2

− εβ(
µ
(
x2 + y2

)
+ z2

)p/2 , (1.6)

where µ is near 1, β is a constant, p is an arbitrary integer and ε is a small parameter.
Originally, in quantum mechanics, the parameter µ of (1.6) is taken as µ 6= 1 and
the problem is called a spatial anisotropic problem.

In this paper, we always assume that µ = 1− ε. As usual, firstly we use a series
of canonical transformations to transform (1.6) into the spatial Delaunay variables
(see [5, 38])

H = − 1

2L2
+ εH1(l, g, L,G,K) +O(ε2), (1.7)

where H1(l, g, L,G,K) is equal to

2

4L
(
L−
√
L2 −G2 cosE

) − β (L(L−√L2 −G2 cosE
))−p

−
2
(
G2 −K2

) ((
L cosE −

√
L2 −G2

)
sin g +G sinE cos g

)2
4G2L

(
L−
√
L2 −G2 cosE

)3
and E = E(l, L,G) is defined by l = E − e sinE, l is the mean anomaly, g is the
argument of the perigee of the unperturbed elliptic orbit measured in the invariant
plane, L is the square root of the semimajor axis of the unperturbed elliptic orbit,
G is the modulus of the total angular momentum, and K is the third component
of the angular momentum.

To state our main results, for any fixed constant h, p, we define two functions as
follows

D1(G;h, p) = −β
√
−2h

p−3
G−p−1

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

×
(

(1− p)(1−G
√
−2h)−m

)(1−G
√
−2h

G
√
−2h

)m−1

and

D2(G;h, p) = −β2p−1
√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(1− p)(1−G
√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1
.

The existence of periodic orbits of Hamiltonian system corresponding to (1.6) is
determined by the functions D1(G;h, p) and D2(G;h, p). Considering the calcula-
tion of abnormal integrals for averaging systems, we divided into two cases, p > 1
and p ≤ 1. Now we state the first result as follows.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 5

Theorem 1.1. Assume that p > 1. For every k0 ∈ [0, 2π), there exists a sufficiently
small positive constant ε0 such that, for each ε ∈ (−ε0, ε0), on every energy level
H = h < 0, the spatial anisotropic Kepler problem with anisotropic perturbations
given by the Hamiltonian (1.7) has

(i) Two 2π-periodic solutions γε
±(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if

|β| < 1

2p+1p(p− 1)|h|p−1
,

such that

lim
ε→0

γε
±(0) =

(
± 1

2
arccos

(
−β2p+1p(p− 1)(−h)p−1

)
, k0,

1√
−2h

,
1√
−2h

, 0
)
.

(ii) Four 2π-periodic solutions γε
(i)(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if

β ∈
(
− 1

2p+1(−h)p−1p(p− 1)
, 0
)
,

and G̃i are two nondegenerate zeros of D1(G;h, p) +
(
2G
√
−h+

√
2
)−2

,
such that

lim
ε→0

γε
(i)(0) =

(
g0, k0,

1√
−2h

, G̃i, 0
)
, g0 = 0 or π, i = 1, 2.

(iii) Four 2π-periodic solutions γε
(i)(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if

β ∈
(

0,
1

2p+1(−h)p−1p(p− 1)

)
,

and Gi are two nondegenerate zeros of D1(G;h, p) −
(
2G
√
−h+

√
2
)−2

,
such that

lim
ε→0

γε
(i)(0) =

(
± π

2
, k0,

1√
−2h

,Gi, 0
)
, i = 1, . . . , 4.

Theorem 1.1 establishes a generic result on the existence of periodic orbits of
the spatial anisotropic Kepler problem with anisotropic perturbations, which gen-
eralizes the results of Lliber and Makhlouf in [29], where the only case p = 2 has
been considered. The proof of Theorem 1.1 combines the Hamiltonian perturbation
theory, Residue theory in complex analysis and the average theory. In the proof of

Theorem 1.1, we have proved the functions D1(G;h, p) ±
(
2G
√
−h+

√
2
)−2

have
at least two different zeros. Owing to lots of calculations, the proof is so long that
we arrange the proof of Theorem 1.1 in Subsection 6.

Similarly, we can obtain the following result for the case p ≤ 1.

Theorem 1.2. Assume that p ≤ 1. For every k0 ∈ [0, 2π), there exists a sufficiently
small positive constant ε0 such that, for each ε ∈ (−ε0, ε0), on every energy level
H = h < 0, the spatial anisotropic Kepler problem with anisotropic perturbations
given by the Hamiltonian (1.7) has

(i) Two 2π-periodic solutions γε
±(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if

|β| < 1

2p+1(p− 1)p|h|p−1
,

such that

lim
ε→0

γε
±(0) =

(
± 1

2
arccos

(
β2p+1(1− p)p(−h)p−1

)
, k0,

1√
−2h

,
1√
−2h

, 0
)
.



6 M. LI, Q. LIU EJDE-2021/63

(ii) Two 2π-periodic solutions γε(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if
p < 0 ,

β ∈
(
−∞,− 1

2p+1(p− 1)p(−h)p−1

)
,

and G̃ is a nondegenerate zero of the D2(G;h, p)+
(
2G
√
−h+

√
2
)−2

, such
that

lim
ε→0

γε
(i)(0) =

(
g0, k0,

1√
−2h

, G̃, 0
)
, g0 = 0 or π, i = 1, 2.

(iii) Two 2π-periodic solutions γε
±(l) = (g(l, ε), k(l, ε), L(l, ε), G(l, ε),K(l, ε)) if

p < 0,

β ∈
( 1

2p+1(p− 1)p(−h)p−1
,+∞

)
,

and G is a nondegenerate zero of the D2(G;h, p)−
(
2G
√
−h+

√
2
)−2

, such
that

lim
ε→0

γε
±(0) =

(
± π

2
, k0,

1√
−2h

,G, 0
)
.

The proof of Theorem 1.2 shows that D2(G;h, p) ±
(
2G
√
−h+

√
2
)−2

have no
zero points for the case of p = 0, 1, while for p < 0, these functions have at least
one zero point. The case p ≤ 0 is corresponding to the anisotropic perturbations
without singularities. We arrange the proof of Theorem 1.2 in Subsection 6.1.

The rest paper is organized as follows. In Section 2, we give some applica-
tions related to spatial anisotropic problem. To demonstrate the applications of
Theorems 1.1 and 1.2, we consider two cases p = 2 and p = −1. The numerical
simulations of periodic orbits are shown in configuration space xyOz. In Section 3,
we introduce two canonical transformations. One is the spherical transformation
Ψ0 : (x, y, z,X, Y, Z)→ (ρ, θ, φ, P,Θ,Φ), and the other one is the spatial Delaunay
transformation Ψ1 : (ρ, θ, φ, P,Θ,Φ) → (l, g, k, L,G,K), which transforms Hamil-
tonian (1.6) into Hamiltonian (1.7) with the Delaunay coordinates. In Section 4, we
obtain the Hamiltonian system corresponding to the Hamiltonian (1.7), and reduce
this system to a two dimensional system on the energy surface Ω, see (4.4). We
emphasize on that we are different from [29] in methodology in this content. In Sec-
tion 5, we use the Residue theorem to obtain the average system of the Hamiltonian
system corresponding to the Hamiltonian (1.7) for p > 1 and p ≤ 1, respectively.
This is the preparation for the next section to use the averaging method to find
the equilibrium points of system. In Section 6, we use averaging theory to prove
Theorem 1.1 and Theorem 1.2, and calculate the equilibrium points of the average
system (5.7) in three cases.

2. Applications related to spatial anisotropic problems

In this section, we give some applications related to spatial anisotropic problem.
The main purpose is to demonstrate the applications of Theorems 1.1 and 1.2. We
shall consider two cases p = 2 and p = −1. The first case is corresponding to
the spatial anisotropic Manev problem. Although this situation has been studied
in [29], only the transformed system is considered. We will return the result of
Theorem 1.1 (see Theorem 2.1 below) to the original system with the Cartesian
coordinates. Moreover, some numerical examples will be exhibited. The second



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 7

case is corresponding to the spatial anisotropic Kepler problem with anisotropic
perturbations.

2.1. The spatial anisotropic Manev problem. When p = 2, we have the spatial
anisotropic Manev problem, which is given by the Hamiltonian

H =
1

2

(
X2 + Y 2 + Z2

)
− 1√

µ
(
x2 + y2

)
+ z2

− εβ

µ
(
x2 + y2

)
+ z2

. (2.1)

Without loss of generality, we consider the periodic orbits of the spatial anisotropic
Manev problem on the energy level H = h = −1. By using Theorem 1.1, we have
the following result.

Theorem 2.1. For every k0 ∈ [0, 2π), there exists a sufficiently small positive
constant ε0 such that, for each ε ∈ (−ε0, ε0), on energy level H = h = −1, the
spatial anisotropic Manev problem given by the Hamiltonian (2.1) has:

(i) Two 2π-periodic solutions ζε
±(t) = (x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε),

Z(t, ε)) if |β| < 1/16, such that ζε
±(0) tends to(1

2
,

1

2
cos g0 sin k0,

1

2
sin g0,

√
2 sin g0 cos k0,

√
2 sin g0 sin k0,−

√
2 cos g0

)
,

as ε→ 0 where

g0 = ±1

2
arccos (−16β) .

(ii) Four 2π-periodic solutions

ζε
i(t) =

(
x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε), Z(t, ε)

)
if β ∈

(
− 1/16, 0

)
, such that ζε

i(0) tends to(
± 1

2

(
1−

√
1− 2G2

0

)
cos k0,±

1

2

(
1−

√
1− 2G2

0

)
sin k0, 0, 0, 0,

∓2|G0|
1−

√
1− 2G2

0

)
,

as ε→ 0 where

G0 =
2
√

2β ±
√
−2β

−1− 4β
.

(iii) Four 2π-periodic solutions

ζε
i(t) =

(
x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε), Z(t, ε)

)
if β ∈ (0, 1/16), such that ζε

i(0) tends to(
0, 0,±1

2

(
1−

√
1− 2G2

0

)
,± 2|G0| cos k0

1−
√

1− 2G2
0

,± 2|G0| sin k0
1−

√
1− 2G2

0

, 0
)
,

as ε→ 0 where

G0 =
2
√

2β ±
√

2β

1− 4β
.

Proof. When |β| < 1/16, from the statement (i) of Theorem 1.1, we obtain

lim
ε→0

γε
±(0) =

(
± 1

2
arccos (−16β) , k0,

1√
2
,

1√
2
, 0
)
,

or equivalently,

γ±0 = (l0, g0, k0, L0, G0,K0) =
(

0,±1

2
arccos (−16β) , k0,

1√
2
,

1√
2
, 0
)
.



8 M. LI, Q. LIU EJDE-2021/63

We insert p = 2 and h = −1 back into the function D1(G;h, p)+
(
2G
√
−h+

√
2
)−2

when β ∈
(
− 1

16 , 0
)
, then obtain

D1(G;h, p) +
(

2G
√
−h+

√
2
)−2

=
G2 + 2β

(
G
√

2 + 1
)2

2G2
(
G
√

2 + 1
)2 . (2.2)

By direct computation, equation (2.2) has two zero points as

G0 =
2
√

2β ±
√
−2β

−1− 4β
,

and these zero points are nondegenerate, i.e.,

d

dG

(
D1(G0) +

(
2G0 +

√
2
)−2 )

6= 0.

From statement (ii) of Theorem 1.1, we obtain

lim
ε→0

γε
i(0) =

(
g0, k0,

1√
2
,

2
√

2β ±
√
−2β

−1− 4β
, 0
)
, g0 = 0 or π, i = 1, . . . , 4,

or equivalent to

γi0 = (l0, g0, k0, L0, G0,K0) =
(

0, g0, k0,
1√
2
,

2
√

2β ±
√
−2β

−1− 4β
, 0
)
,

g0 = 0 or π, i = 1, . . . , 4.

We insert p = 2 and h = −1 back into D1(G;h, p) −
(
2G
√
−h+

√
2
)−2

when
β ∈ (0, 1/16), then obtain

D1(G;h, p)−
(

2G
√
−h+

√
2
)−2

=
−G2 + 2β

(
G
√

2 + 1
)2

2G2
(
G
√

2 + 1
)2 . (2.3)

The zero points of (2.3) are

G0 =
2
√

2β ±
√

2β

1− 4β
.

We derive D1(G) −
(
2G+

√
2
)−2

and replace G in the derived equation with G0,
obtain

d

dG

(
D1(G0)−

(
2G0 +

√
2
)−2 )

6= 0.

Therefore, G0 are the nondegenerate zero of the D1(G)−
(
2G+

√
2
)−2

. From the
statement (iii) of Theorem 1.1, we obtain

lim
ε→0

γε
i(0) =

(
± π

2
, k0,

1√
2
,

2
√

2β ±
√

2β

1− 4β
, 0
)
, i = 1, . . . , 4,

or equivalent to

γi0 = (l0, g0, k0, L0, G0,K0) =
(

0,±π
2
, k0,

1√
2
,

2
√

2β ±
√

2β

1− 4β
, 0
)
, i = 1, . . . , 4.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 9

To study periodic solutions of the system corresponding to Hamiltonian (2.1), we
need two to transform the Delaunay coordinates (l, g, k, L,G,K) into the Cartesian
coordinates (x, y, z,X, Y, Z). Firstly, inserting

(l0, k0, L0,K0) =
(

0, k0,
1√
2
, 0
)
, k0 ∈ [0, 2π)

into the Delaunay transformation Ψ1 (see (3.3)), and transforming the Delaunay
coordinates (l, g, k, L,G,K) into the spatial spherical coordinates (ρ, θ, φ, P,Θ,Φ),
we have

ρ =
1

2

(
1−

√
1− 2G2

)
, θ = k0, φ =

π

2
− g, P = 0,Θ = 0,Φ = |G|. (2.4)

Then we put (2.4) into the spherical transformation Ψ0 (see (3.1)) and obtain

x =
1

2

(
1−

√
1− 2G2

)
cos g cos k0,

y =
1

2

(
1−

√
1− 2G2

)
cos g sin k0,

z =
1

2

(
1−

√
1− 2G2

)
sin g,

cos g (X cos k0 + Y sin k0) + Z sin g = 0,

1

2

(
1−

√
1− 2G2

)
cos g (−X sin k0 + Y cos k0) = 0,

1

2

(
1−

√
1− 2G2

)
[sin g (X cos k0 + Y sin k0)− Z cos g] = |G|.

(2.5)

Considering the statement (i) of Theorem 1.1, we insert

(g0, G0) =
(
± 1

2
arccos (−16β) ,

1√
2

)
into (2.5) to obtain

ζ±0 = (x0, y0, z0, X0, Y0, Z0)

=
(1

2
cos g0 cos k0,

1

2
cos g0 sin k0,

1

2
sin g0,

√
2 sin g0 cos k0,

√
2 sin g0 sin k0,−

√
2 cos g0

)
.

Similarly, considering the statement (ii) of Theorem 1.1, we take

(g0, G0) =
(

0 or π,
2
√

2β ±
√
−2β

−1− 4β

)
into (2.5) to obtain

ζi0 =
(
± 1

2

(
1−

√
1− 2G2

0

)
cos k0,±

1

2

(
1−

√
1− 2G2

0

)
sin k0,

0, 0, 0,∓ 2|G0|
1−

√
1− 2G2

0

)
, i = 1, . . . , 4.

Considering the statement (iii) of Theorem 1.1, we take

(g0, G0) =
(
± π

2
,

2
√

2β ±
√

2β

1− 4β

)



10 M. LI, Q. LIU EJDE-2021/63

into (2.5) to obtain

(x0, y0, z0, Z0) =
(

0, 0,±1

2

(
1−

√
1− 2G2

0

)
, 0
)
, (2.6)

±1

2

(
1−

√
1− 2G2

0

)
(X0 cos k0 + Y0 sin k0) = |G0|. (2.7)

Inserting (2.6) and H = h = −1 back into Hamiltonian (2.1), we obtain

X2
0 + Y 2

0 = −2 +
4(1−

√
1− 2G2

0) + 8εβ

(1−
√

1− 2G2
0)2

. (2.8)

When ε→ 0, combining (2.7) with (2.8), we obtain

(X0, Y0) =
(
± 2|G0| cos k0

1−
√

1− 2G2
0

,± 2|G0| sin k0
1−

√
1− 2G2

0

)
.

Therefore,

ζi0 =
(

0, 0,±1

2

(
1−

√
1− 2G2

0

)
,± 2|G0| cos k0

1−
√

1− 2G2
,± 2|G0| sin k0

1−
√

1− 2G2
0

, 0
)
,

for i = 1, . . . , 4. Thus the proof is complete. �

From Theorem 2.1, when β ∈ (0, 1/16) or β ∈ (−1/16, 0), for sufficiently small
|ε| � 1, the spatial anisotropic Manev problem corresponding to Hamiltonian (2.1)
have at least six kinds of periodic solutions, which constitute a continuum of peri-
odic solutions by continuously changing k0 on the interval [0, 2π).

To demonstrate Theorem 2.1, we take the parameter β = 1/32 and k0 = π/3.
According to Theorem 2.1, there exists six periodic orbits ζiε(t) such that

lim
ε→0

ζε
i(0) = ζi(0)

=



(
1
8 ,
√
3
8 ,
√
3
4 ,
√
6
4 ,

3
√
2

4 ,− 1√
2

)
for i = 1,(

1
8 ,
√
3
8 ,−

√
3
4 ,−

√
6
4 ,−

3
√
2

4 ,− 1√
2

)
for i = 2,(

0, 0, 12 −
√

10−
√
2

7 ,−
√
2
2 +

√
4− 2

√
2 + 2,−

√
6
2 +

√
12− 6

√
2 + 2

√
3, 0
)
,(

0, 0, 12 −
√

10+
√
2

7 ,
√
2
2 +

√
4 + 2

√
2 + 2,

√
6
2 +

√
12 + 6

√
2 + 2

√
3, 0
)
,(

0, 0,

√
10−
√
2

7 − 1
2 ,
√
2
2 −

√
4− 2

√
2− 2,

√
6
2 −

√
12− 6

√
2− 2

√
3, 0
)
,(

0, 0,

√
10+
√
2

7 − 1
2 ,−

√
2
2 −

√
4 + 2

√
2− 2,−

√
6
2 −

√
12 + 6

√
2− 2

√
3, 0
)
.

For numerical demonstration of periodic orbits, we take ε = 0.002 and µ = 0.998.
Since for sufficiently small |ε| � 1, the initial values of periodic solutions ζiε(t), i =
1, 2, . . . , 6 are near ζi(0), we plot the curves of the solutions starting at ζi(0) on the
time interval t ∈ (990, 1000), see Figure 1. As is shown in Figure 1, the solutions
starting at ζi(0) are periodic. Since the flows of autonomous Hamiltonian system
have the semigroup property, two different periodic solutions may be have the same
orbits. In order to test different periodic orbits, we plot the periodic orbits ζiε(t)
(i = 1, 2, . . . , 6) on the projective configuration space xOyz, see Figure 2(a). As can
be seen from Figure 2(a), there are five differential periodic orbits, since periodic



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 11

Figure 1. Periodic orbits of Hamiltonian (2.1) with the parame-
ters ε = 0.002, β = 1/32 and k0 = π/3. In the i-th (i = 1, 2, . . . , 6)
line, the curves represent periodic orbit ζiε(t) corresponding to the
initial value ζi(0), respectively. On the left, the pictures denote by
the curves x(t), y(t), z(t), while on the right, the pictures denote
by the curves X(t), Y (t), Z(t).

solutions ζ1ε (t) and ζ2ε (t) have the same periodic orbit. Moreover, we plot the orbits
on the projective configuration space xOyz starting from the initial values

ζ1(0) =
(1

4
cos k0,

1

4
sin k0,

√
3

4
,

√
3

2
cos k0,

√
3

2
sin k0,−

1√
2

)
by varying k0 = iπ/6, i = 1, 2, . . . , 6, see Figure 2(b). As is shown from Figure 2(b),
if we continuously change k0 on the interval (0, 2π), all of periodic orbits are full of
an ellipsoid.

2.2. The spatial anisotropic Kepler problem with perturbations. In this
subsection, we consider the spatial anisotropic Kepler problem with anisotropic
perturbations. We only consider p = −1 as an example for the case p ≤ 1. In
this case, the integral computation of the perturbations does not need the Residue



12 M. LI, Q. LIU EJDE-2021/63

k0 = π/3 and (xiε(t), y
i
ε(t), z

i
ε(t)) k0 = iπ/6 (i = 1, 2, . . . , 6)

Figure 2. Periodic orbits on the configuration space xOyz,
where we choose the parameters as ε = 0.002 and β = 1/32. We
plot the periodic orbits ζiε(t) (i = 1, 2, . . . , 6) on the projective
configuration space xOyz on the left, while on the right we
plot the periodic orbits starting from the initial values ζ1(0) =

(cos k0/4, sin k0/4,
√

3/4,
√

3/2 cos k0,
√

3/2 sin k0,−1/
√

2) by
varying k0 = iπ/6, i = 1, 2, . . . , 6.

theorem. The Hamiltonian is

H =
1

2

(
X2 + Y 2 + Z2

)
− 1√

µ
(
x2 + y2

)
+ z2

− εβ
√
µ
(
x2 + y2

)
+ z2. (2.9)

Without loss of generality, we consider the periodic orbits of the spatial anisotropic
Kepler problem with anisotropic perturbations on the energy level H = h = −1.
By using Theorem 1.2, we have the following result.

Theorem 2.2. For every k0 ∈ [0, 2π), there exists a sufficiently small positive con-
stant ε0 such that, for each ε ∈ (−ε0, ε0), on energy level H = h = −1, the spatial
anisotropic Kepler problem with anisotropic perturbations defined by the Hamilton-
ian (2.9) has:

(i) Two 2π-periodic solutions

ζε
±(t) = (x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε), Z(t, ε))

such that if |β| < 1/2, then ζε
±(0) tends to(1

2
cos g0 cos k0,

1

2
cos g0 sin k0,

1

2
sin g0,

√
2 sin g0 cos k0,

√
2 sin g0 sin k0,−

√
2 cos g0

)
,

as ε→ 0, where

g0 = ±1

2
arccos (−2β) .

(ii) Two 2π-periodic solutions

ζε
±(t) = (x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε), Z(t, ε))



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 13

such that if β ∈ (−∞,−1/2), then ζε
±(0) tends to(

± 1

2

(
1−

√
1− 2G2

0

)
cos k0,±

1

2

(
1−

√
1− 2G2

0

)
sin k0, 0, 0, 0,∓

2|G0|
1−

√
1− 2G2

0

)
,

as ε→ 0, where

G0 = −
3

√
−β2

(
β + 3

(√
81− 6β − 9

))
3
√

2β
+

1

3
√

2
3

√
β

β + 3
(√

81− 6β − 9
) − √2

3
.

(iii) Two 2π-periodic solutions

ζε
±(t) = (x(t, ε), y(t, ε), z(t, ε), X(t, ε), Y (t, ε), Z(t, ε))

such that if β ∈ (1/2,+∞), then ζε
±(0) tends to(

0, 0,±1

2

(
1−

√
1− 2G2

0

)
,± 2|G0| cos k0

1−
√

1− 2G2
0

,± 2|G0| sin k0
1−

√
1− 2G2

0

, 0
)
,

as ε→ 0, where

G0 =

(
β −

(
3
√

3
√
β4(2β + 27) + (β + 27)β2

)1/3)2
3
√

2β
(
3
√

3
√
β4(2β + 27) + (β + 27)β2

)1/3 .

Proof. By using Theorem 2.2, in each case we can obtain two periodic solutions
of Hamiltonian (2.9) such that limε→0 γε

±(0) = (g0, k0, L0, G0,K0). Then we use
the Delaunay transformation Ψ1 (see (3.3)) to change the Delaunay coordinates
(0, g0, k0, L0, G0,K0) into the spatial spherical coordinates (ρ0, θ0, φ0, P0,Θ0,Φ0).
Subsequently, we use the spherical transformation Ψ0 (see (3.1)) to transform the
spatial spherical coordinates into the Cartesian coordinates (x0, y0, z0, X0, Y0, Z0).
The rest proof of Theorem 2.2 is similar to the one of Theorem 2.1. To avoid
repetition, we do not do it. �

Figure 3. Periodic orbits of Hamiltonian (2.9) with the param-
eters ε = 0.002, β = −1 and k0 = π/3. In the i-th (i = 1, 2)
line, the curves represent periodic orbit ζiε(t) corresponding to the
initial value ζi(0), respectively. On the left, the pictures denote by
the curves x(t), y(t), z(t), while on the right, the pictures denote
by the curves X(t), Y (t), Z(t).

From Theorem 2.2, when β ∈ (−∞,−1/2), β ∈ (−1/2, 0) or β ∈ (0, 1/2),
β ∈ (1/2,+∞), for sufficiently small |ε| � 1, the spatial anisotropic Kepler problem
with perturbations corresponding to Hamiltonian (2.9) have at least two kinds of



14 M. LI, Q. LIU EJDE-2021/63

k0 = π/3 and (xiε(t), y
i
ε(t), z

i
ε(t)) k0 = iπ/6 (i = 1, 2, . . . , 6)

Figure 4. Periodic orbits on the configuration space xOyz,
where we choose the parameters as ε = 0.002 and β = −1.
We plot the periodic orbits ζiε(t) (i = 1, 2) on the projective
configuration space xOyz on the left, while on the right we
plot the periodic orbits starting from the initial values ζ1(0) =
(0.140795 cos k0, 0.140795 sin k0, 0, 0,−3.49357) by varying k0 =
iπ/6, i = 1, 2, . . . , 6.

periodic solutions, which constitute a continuum of periodic solutions by continu-
ously changing k0 on the interval [0, 2π).

To demonstrate Theorem 2.2, we take the parameter β = −1 and k0 = π/3.
According to Theorem 2.2, there exists two periodic orbits ζiε(t), i = 1, 2 such that

lim
ε→0

ζε
i(0) = ζi(0) =

{
(0.070398, 0.121932, 0, 0, 0,−3.493567)

(−0.070398,−0.121932, 0, 0, 0, 3.493567) .

For numerical demonstration of periodic orbits, we take ε = 0.002 and µ =
0.998. Since for sufficiently small |ε| � 1, the initial values of periodic solutions
ζiε(t), i = 1, 2 are near ζi(0), we plot the curves of the solutions starting at ζi(0)
on the time interval t ∈ (990, 1000), see Figure 3. As is shown in Figure 3, the
solutions starting at ζi(0) are periodic. In order to test different periodic orbits, we
plot the periodic orbits ζiε(t) (i = 1, 2) on the projective configuration space xOyz,
see Figure 4(a). As can be seen from Figure 4(a), there are two differential periodic
orbits. Moreover, we plot the orbits on the projective configuration space xOyz
starting from the initial values

ζ1(0) = (0.140795 cos k0, 0.140795 sin k0, 0, 0, 0,−3.493567)

by varying k0 = iπ/6, i = 1, 2, . . . , 6, see Figure 4(b). As shown IN Figure 4(b), if
we continuously change k0 on the interval [0, 2π), all of periodic orbits are full of
an ellipsoid.

Similarly, to demonstrate Theorem 2.2, we take the parameter β = 1/32 and
k0 = π/3. According to Theorem 2.2, there exists two periodic orbits ζiε(t), i = 1, 2
such that

lim
ε→0

ζε
i(0) = ζi(0) =


( √

15
16
√
2
, 3
√
5

16
√
2
,
√
17

8
√
2
,
√
17
8 ,

√
51
8 ,−

√
15
4

)
for i = 1,( √

15
16
√
2
, 3
√
5

16
√
2
,−
√
17

8
√
2
,−
√
17
8 ,−

√
51
8 ,−

√
15
4

)
for i = 2.

For numerical examples of periodic orbits, we take ε = 0.002 and µ = 0.998.
Since for sufficiently small |ε| � 1, the initial values of periodic solutions ζiε(t), i =



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 15

Figure 5. Periodic orbits of Hamiltonian (2.9) with the param-
eters ε = 0.002, β = 1/32 and k0 = π/3. In the i-th (i = 1, 2)
line, the curves represent periodic orbit ζiε(t) corresponding to the
initial value ζi(0), respectively. On the left, the pictures denote by
the curves x(t), y(t), z(t), while on the right, the pictures denote
by the curves X(t), Y (t), Z(t).

k0 = π/3 and (xiε(t), y
i
ε(t), z

i
ε(t)) k0 = ipi/6 (i = 1, 2, . . . , 6)

Figure 6. Periodic orbits on the configuration space xOyz,
where we choose the parameters as ε = 0.002 and β = 1/32.
We plot the periodic orbits ζiε(t) (i = 1, 2) on the projec-
tive configuration space xOyz on the left, while on the right
we plot the periodic orbits starting from the initial values
ζ1(0) = (

√
15 cos k0/8

√
2,
√

15 sin k0/8
√

2,
√

17/8
√

2,
√

17 cos k0/4,√
17 sin k0/4,−

√
15/4) by varying k0 = iπ/6, i = 1, 2, . . . , 6.

1, 2 are near ζi(0), we plot the curves of the solutions starting at ζi(0) on the time
interval t ∈ (990, 1000), see Figure 5. As is shown in Figure 5, the solutions starting
at ζi(0) are periodic. Since the flows of autonomous Hamiltonian system have the
semigroup property, two different periodic solutions may be have the same orbits.
In order to test different periodic orbits, we plot the periodic orbits ζiε(t) (i = 1, 2)
on the projective configuration space xOyz, see Figure 6(a). As can be seen from
Figure 6(a), there is only one periodic orbit, since periodic solutions ζ1ε (t) and
ζ2ε (t) have the same periodic orbit. Moreover, we plot the orbits on the projective
configuration space xOyz starting from the initial values

ζ1(0) =
(√15

8
√

2
cos k0,

√
15

8
√

2
sin k0,

√
17

8
√

2
,

√
17

4
cos k0,

√
17

4
sin k0,−

√
15

4

)



16 M. LI, Q. LIU EJDE-2021/63

by varying k0 = iπ/6, i = 1, 2, . . . , 6, see Figure 6(b). As is shown from Figure 6(b),
if we continuously change k0 on the interval (0, 2π), all of periodic orbits are full of
an ellipsoid.

3. Canonical transformations

In this section, we introduce two canonical transformations, which will be used
to transform Hamiltonian (1.6) with the Cartesian coordinates (x, y, z,X, Y, Z) into
Hamiltonian (1.7) with the Delaunay coordinates (l, g, k, L,G,K). The content of
this part is known and we can refer to some classical books (for example, see [5, 38])
for details.

3.1. Spherical coordinates. Firstly, we shall introduce a symplectic transforma-
tions

Ψ0 : (x, y, z,X, Y, Z)→ (ρ, θ, φ, P,Θ,Φ)

by Ψ0:

x = ρ cos θ sinφ, y = ρ sin θ sinφ, z = ρ cosφ,

X =
cos θ(Pρ sinφ+ Φ cosφ)−Θ sin θ cscφ

ρ
,

Y =
Θ cos θ cscφ+ sin θ(Pρ sinφ+ Φ cosφ)

ρ
,

Z = P cosφ− Φ sinφ

ρ
,

(3.1)

which is generated by the Mathieu generating function

W0(ρ, θ,X, Y, Z) = ρX cos θ sinφ+ ρY sin θ sinφ+ ρZ cosφ.

Under the symplectic transformations Ψ0, the Hamiltonian (1.6) becomes

H =
1

2

(Θ2 csc2(φ)

ρ2
− 2

√
2√

ρ2(1 + µ+ (1− µ) cos(2φ))
+ P 2 +

Φ2

ρ2

− ε2
p
2+1β

(
ρ2
(

1 + µ+ (1− µ) cos(2φ)
))−p/2 )

.

(3.2)

3.2. Spatial Delaunay elements. To change from the spherical coordinates
(ρ, θ, φ, P,Θ,Φ) to the Delaunay elements (l, g, k, L,G,K) with the first three vari-
ables defined mod 2π, we consider the generating function

W (ρ, θ, φ, L,G,K) = θK +

∫ φ

π/2

(
G2− K2

sin2 ξ

)1/2
dξ+

∫ ρ

ρ0

(
− G2

ξ2
− 1

L2
+

2

ξ

)1/2
dξ,

where ρ0 = L2(1−
√

1−G2/L2).



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 17

Then the change of coordinates is Ψ1:

P =
(
− G2

ρ2
− 1

L2
+

2

ρ

)1/2
,

Θ = K,

Φ =
(
G2 − K2

sin2 φ

)1/2
,

l =
1

L3

∫ ρ

ρ0

(
− G2

ξ2
− 1

L2
+

2

ξ

)−1/2
dξ,

g = −
∫ φ

π/2

(
G2 − K2

sin2 ξ

)−1/2
Gdξ −

∫ ρ

ρ0

(
− G2

ξ2
− 1

L2
+

2

ξ

)−1/2(G
ξ2

)
dξ,

k = θ −
∫ φ

π/2

(
G2 − K2

sin2 ξ

)−1/2( K2

sin2 ξ

)
dξ.

(3.3)

The first integral in (3.3) can be obtained explicitly by

l =
1

L3

∫ ρ

ρ0

(
− G2

ξ2
− 1

L2
+

2

ξ

)−1/2
dξ = E − e sinE, (3.4)

where E is determined by ρ = a(1 − e cosE) with a = L2 and e =
√

1−G2/L2.
Thus we can see that the variable ρ = ρ(l, L,G) depends on l, L,G. We use −σ, f
to denote the second integral and the third integral, respectively, where cosφ =
sinσ sin i = sin(f + g) sin i with

cos i = K/G, ρ cos f = a(cosE − e), ρ sin f =
a

L
G sinE.

From the equalities above, we know that the variable φ = φ(l, g, L,G,K) depends
on l, g, L,G,K. The last equation of (3.3) implies that sin(k−θ) = (cotφ)/γ, γ2 =
(G2 −K2)/K2.

Under the symplectic transformations Ψ1, the Hamiltonian (3.2) becomes

H = − 1

2L2
+

1

ρ
−

√
2

ρ
√

(1 + µ+ (1− µ) cos(2φ))

− ε2
p
2 β
[
ρ2
(

1 + µ+ (1− µ) cos(2φ)
)]−p/2

,

(3.5)

where ρ = ρ(l, L,G), φ = φ(l, g, L,G,K). We emphasize that the Hamiltonian (3.5)
is independent of the variable k, which implies that K(t) is always constant.

Let µ = 1− ε. We expand the Hamiltonian function (3.5) at ε = 0 into

H =− 1

2L2
− ε
(
βρ−p +

1− cos(2φ)

4ρ

)
− 1

8
ε2
(

2βpρ−p
(
1− cos(2φ)

)
+

3
(
3− 4 cos(2φ) + cos(4φ)

)
8ρ

)
+O(ε3).

(3.6)

Since

ρ2 cos2 φ = ρ2 sin2(f + g) sin2 i = (ρ sin f cos g + ρ cos f sin g)
2

(1− cos2 i),



18 M. LI, Q. LIU EJDE-2021/63

by using the relationships between various variables, we know that

ρ2 cos(2φ) =2 (ρ sin f cos g + ρ cos f sin g)
2

(1− cos2 i)− ρ2

=2
( a
L
G sinE cos g + a(cosE − e) sin g

)2 (
1− K2

G2

)
− ρ2.

(3.7)

Substituting (3.7) and ρ = a(1− e cosE) into (3.6), we obtain

H =− 1

2L2
+ ε

(
1

2L
(
L−
√
L2 −G2 cosE

) − β (L(L−√L2 −G2 cosE
))−p

−
(G2 −K2)

(
(L cosE −

√
L2 −G2) sin g +G sinE cos g

)2
2G2L(L−

√
L2 −G2 cosE)3

)
+O(ε2)

:=H0(L) + εH1(l, g, L,G,K) +O(ε2),

(3.8)
where E = E(l, L,G) is defined by l = E − e sinE.

4. Reduction of system on energy surface

The Hamiltonian system corresponding to the Hamiltonian (3.8) is given by

dl

dt
=
∂H

∂L
=

1

L3
+O(ε),

dg

dt
=
∂H

∂G
= ε

∂H1

∂G
(l, g, L,G,K) +O(ε2),

dk

dt
=
∂H

∂K
= ε

∂H1

∂K
(l, g, L,G,K) +O(ε2),

dL

dt
= −∂H

∂l
= −ε∂H1

∂l
(l, g, L,G,K) +O(ε2),

dG

dt
= −∂H

∂g
= −ε∂H1

∂g
(l, g, L,G,K) +O(ε2),

dK

dt
= −∂H

∂k
= 0,

(4.1)

where the last equation is owing to that H is independent of k, see (3.5).
For sufficiently small |ε| � 1 and L > 0, we know that dl/dt = 1/L3 +O(ε) > 0.

Then we regard the angular variable l as a new time variable. From (4.1), we can
obtain a new system of differential equations

dg

dl
= εL3 ∂H1

∂G
(l, g, L,G,K) +O(ε2),

dk

dl
= εL3 ∂H1

∂K
(l, g, L,G,K) +O(ε2),

dL

dl
= −εL3 ∂H1

∂l
(l, g, L,G,K) +O(ε2),

dG

dl
= −εL3 ∂H1

∂g
(l, g, L,G,K) +O(ε2),

dK

dl
= −∂H

∂k
= 0.

(4.2)

Since l is an angular variable defined mod 2π, H,H1 and system (4.2) are also
2π-periodic in l. For any h < 0, we define the energy surface by

Ω =
{

(l, g, k, L,G,K) ∈ R6 : H(l, g, k, L,G,K) = h
}
. (4.3)



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 19

For sufficiently small |ε| � 1, we obtain that L = 1/
√
−2h+O(ε). By using the first

order expansion of L, we can discard the third equation in system (4.2). Moreover,
from the last equation of (4.2) we have K(l) ≡ K0, where K0 is a constant to be
determined later to guarantee that k(l) is 2π-periodic, since the right side of the
second equation of (4.2) dose not depend on k. Therefore, to prove the existence
of periodic solutions of (4.2), we only need to consider

dg

dl
= εL3 ∂H1

∂G
(l, g, L,G,K) +O(ε2),

dG

dl
= −εL3 ∂H1

∂g
(l, g, L,G,K) +O(ε2),

(4.4)

where K = K0 and L = 1/
√
−2h are considered as constants in system (4.4) at

this stage, and the high order part O(ε2) depends on l, g, h,G,K.

5. Standard form of averaging

To investigate the existence of 2π-periodic solutions of system (4.4), we perform
the method of averaging. The averaged system of (4.4) is

dg

dl
= ε

L3

2π

∫ 2π

0

∂H1

∂G
(l, g, L,G,K)dl,

dG

dl
= −εL

3

2π

∫ 2π

0

∂H1

∂g
(l, g, L,G,K)dl.

(5.1)

Since H1, ∂H1/∂G and ∂H1/∂g are continuous, system (5.1) is equivalent to

dg

dl
= ε

L3

2π

∂

∂G

(∫ 2π

0

H1(l, g, L,G,K)dl
)
,

dG

dl
= −εL

3

2π

∂

∂g

(∫ 2π

0

H1(l, g, L,G,K)dl
)
.

(5.2)

To obtain an explicit expression of system (5.2), we have to calculate the definite
integral above. It can be seen from (3.4) that, E varies form 0 to 2π when the
integral variable l varies form 0 to 2π. By the integral method of substitution, we
have

1

2π

∫ 2π

0

H1(l, g, L,G,K)dl =
1

2π

∫ 2π

0

(
H1(l, g, L,G,K)

dl

dE

)
dE. (5.3)

Recall that ρ = L2(1− e cosE), and by (3.4) we have

dl

dE
= 1− e cosE =

ρ

L2
.

For convenience of calculations, we divide H1 · dl/dE into two parts:

H1(l, g, L,G,K)
dl

dE
=

ρ

L2

(
− βρ−p − sin2(φ)

2ρ

)
=
(
− βρ1−p

L2

)
+
(
− sin2(φ)

2L2

)
:= H10(E, g, L,G,K) +H11(E, g, L,G,K).

Since

ρ2 cos2 φ = ρ2 sin2(f + g) sin2 i = (ρ sin f cos g + ρ cos f sin g)
2

(1− cos2 i),



20 M. LI, Q. LIU EJDE-2021/63

by using the relationships between various variables, we know that

ρ2 sin2 φ =ρ2 − (ρ sin f cos g + ρ cos f sin g)
2

(1− cos2 i)

=ρ2 −
( a
L
G sinE cos g + a(cosE − e) sin g

)2 (
1− K2

G2

)
.

(5.4)

Substituting (5.4), a = L2 and ρ = a(1− e cosE) into H10 and H11, we obtain

H10(E, g, L,G,K) =− βρ1−p

L2
= −

β
(
L2(1− e cosE)

)1−p
L2

=− βL−2p

(1− e cosE)p−1

and

H11(E, g, L,G,K)

= − sin2 φ

2L2
=

(sin f cos g + cos f sin g)2G
2−K2

G2 − 1

2L2

=
1

2L2

( (G−K)(G+K)(L sin g(cosE − e) +G sinE cos g)2

G2L2(e cosE − 1)2
− 1
)

= −
sin2E sin2 g

(
G2 −K2

)
4L4(e cosE − 1)2

+
sin2E cos2 g

(
G2 −K2

)
4L4(e cosE − 1)2

+
sinE sin g cos g

(
G2 −K2

)
(cosE − e)

GL3(e cosE − 1)2

−
cos2 g

(
G2 −K2

)
(cosE − e)2

4G2L2(e cosE − 1)2
+

sin2 g
(
G2 −K2

)
(cosE − e)2

4G2L2(e cosE − 1)2

− G2 +K2

4G2L2

:= P1 + P2 + P3 + P4 + P5 + P6.

5.1. Calculation of integral for H10. To calculate the definite integral (5.3), we
calculate the integral for H10

H10 :=
1

2π

∫ 2π

0

H10(E, g, L,G,K) dE = − 1

2π
βL−2p

∫ 2π

0

1

(1− e cosE)p−1
dE.

In the following, by using Residue theory of complex functions we calculate the
definite integral ∫ 2π

0

1

(1− e cosE)p−1
dE

with two different cases for the range of the integer p.

Case I: p > 1. Let z = eiE , then cosE = (z + z−1)/2. When E experiences the
interval [0, 2π], z goes through the unit circle |z| = 1 one round in the positive
direction. Then it follows that∫ 2π

0

1

(1− e cosE)p−1
dE =

1

i

∫
|z|=1

2p−1zp−2

(−ez2 − e+ 2z)
p−1 dz :=

1

i

∫
|z|=1

f(z)dz.

The function

f(z) :=
2p−1zp−2

(−ez2 − e+ 2z)
p−1 =

2p−1zp−2

(−e)p−1
(
z − 1−

√
1−e2
e

)p−1(
z −

√
1−e2+1
e

)p−1



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 21

has only a singular point z0 = (1−
√

1− e2)/e in the open unit disc, which is a pole
of order p − 1. At the same time, the function f is continuous on the unit circle
|z| = 1. Let

g(z) =
2p−1zp−2

(−e)p−1
(
z −

√
1−e2+1
e

)p−1 .
Using the Leibniz expansive formula, we have

dp−2

dzp−2
g(z)

∣∣∣
z=z0

=
2p−1

(−e)p−1
p−2∑
m=0

Cmp−2
(
zp−2

)(−m+p−2)
((
z −
√

1− e2 + 1

e

)1−p)(m)∣∣∣
z=z0

=

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!

(
1−

√
1− e2

)m(√
1− e2

)−m−p+1

=

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!

(
1− G

L

)m(G
L

)−m−p+1
,

where u(m)(z) denotes the m-th order derivatives of u with respect to z, and the
last equality is obtained by using e2 = 1 − G2/L2. By the Residue theorem, we
have∫ 2π

0

1

(1− e cosE)p−1
dE =

1

i

∫
|z|=1

2p−1zp−2

(−ez2 − e+ 2z)
p−1 dz

= 2π Res
z=

1−
√

1−e2

e

f(z)

= 2 π
1

(p− 2)!

dp−2

dzp−2
g(z)

∣∣∣∣
z=

1−
√

1−e2

e

= 2π

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

(
1− G

L

)m(G
L

)−m−p+1

.

Therefore,

H10 = −βL−2p
p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

(
1− G

L

)m(G
L

)−m−p+1

= −βL−p−1
p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!
(L−G)

m
G−m−p+1.

(5.5)

Case II: p ≤ 1. In this case,∫ 2π

0

1

(1− e cosE)p−1
dE =

1

i

∫
|z|=1

(
−ez2 − e+ 2z

)1−p
21−pz2−p

dz :=
1

i

∫
|z|=1

f(z)dz.

We know the singular point z = 0 is a pole of order 2−p. With the same argument
above, we obtain

1

i

∫
|z|=1

(
− ez2 − e+ 2z

)1−p
21−pz2−p

dz



22 M. LI, Q. LIU EJDE-2021/63

=
1

i

∫
|z|=1

(−e)1−p
(
z − 1−

√
1−e2
e

)1−p(
z −

√
1−e2+1
e

)1−p
21−pz2−p

dz

= 2πResz=0 f(z)

= 2π
(−e)1−p

21−p

1−p∑
m=0

Cm1−p

((
z − 1−

√
1− e2
e

)1−p)(1−m−p)
×
((
z −
√

1− e2 + 1

e

)1−p)(m)∣∣∣
z=0

= 2π2p−1
1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

(
1−

√
1− e2

)m(√
1− e2 + 1

)−m−p+1

= 2π2p−1
1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

(
1− G

L

)m(G
L

+ 1
)−m−p+1

.

Therefore,

H10

= −βL−2p2p−1
1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

(
1− G

L

)m(G
L

+ 1
)−m−p+1

= −βL−1−p2p−1
1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!
(L−G)

m
(L+G)

−m−p+1
.

(5.6)

5.2. Calculation of integral for H11. With the same argument for calculating
the integral of H01, we also use the Residue theorem to calculate the integral of
H11. To calculate the integrals of the first two terms P1, P2 of H11, we consider the
following definite integral∫ 2π

0

sin2E

(e cosE − 1)2
dE =

1

i

∫
|z|=1

−
(
z2 − 1

)2
z (ez2 + e− 2z)

2 dz :=
1

i

∫
|z|=1

f1(z)dz,

by setting z = eiE . The function f1(z) has two singular points z1 = 0 and z2 =

(1−
√

1− e2)/e in the open unit disc. Moreover, z1 is a pole of order one while z2
is a pole of order two. Then by using the Residue theorem we have∫ 2π

0

sin2E

(e cosE − 1)2
dE = 2π

2∑
k=1

Res
z=zk

f(z)

= 2π
(
− 1

e2
+

−e2 − 2
√

1− e2 + 2

e2
√

1− e2
(√

1− e2 − 1
)2)

= 2π
1−
√

1− e2

e2
√

1− e2

=
2πL2

G(G+ L)
.

Therefore,

P 1 :=
1

2π

∫ 2π

0

−
sin2E sin2 g

(
G2 −K2

)
4L4(e cosE − 1)2

dE = −
sin2 g

(
G2 −K2

)
4GL2(G+ L)

,



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 23

P 2 :=
1

2π

∫ 2π

0

sin2E cos2 g
(
G2 −K2

)
4L4(e cosE − 1)2

dE =
cos2 g

(
G2 −K2

)
4GL2(G+ L)

.

With the same argument, to calculate the integrals of the first two terms P4, P5;
we consider the definite integral

∫ 2π

0

(cosE − e)2

(e cosE − 1)2
dE =

1

i

∫
|z|=1

(
−2ez + z2 + 1

)2
(ez2 + e− 2z)

2 dz :=
1

i

∫
|z|=1

f2(z)dz.

The function f2(z) has a pole of order one z1 = 0 and a pole of order two z2 =

(1−
√

1− e2)/e. Then by using the Residue theorem we have

∫ 2π

0

(cosE − e)2

(e cosE − 1)2
dE = 2π

2∑
k=1

Res
z=zk

f(z)

= 2π
( 1

e2
−
(
e2 − 1

) (
e2 + 2

√
1− e2 − 2

)
e2
√

1− e2
(√

1− e2 − 1
)2 )

= 2π
1√

1− e2 + 1

=
2πL

G+ L
.

Therefore,

P 4 :=
1

2π

∫ 2π

0

−
cos2 g

(
G2 −K2

)
(cosE − e)2

4G2L2(e cosE − 1)2
dE = −

cos2 g
(
G2 −K2

)
4G2L(G+ L)

,

P 5 :=
1

2π

∫ 2π

0

sin2 g
(
G2 −K2

)
(cosE − e)2

4G2L2(e cosE − 1)2
=

sin2 g
(
G2 −K2

)
4G2L(G+ L)

.

Since the term P3 is an odd function with respect to E, we know the mean value
of P3 along a period 2π is zero, that is, P 3 = 0. Moreover, the term P6 does not
depend on E, then we have the mean value of P6 P 6 = P6.

Finally, we conclude that

H11 : =
1

2π

∫ 2π

0

H11(E, g, L,G,K) dE

= P 1 + P 2 + P 3 + P 4 + P 5 + P 6

= −
sin2 g

(
G2 −K2

)
4GL2(G+ L)

+
cos2 g

(
G2 −K2

)
4GL2(G+ L)

−
cos2 g

(
G2 −K2

)
4G2L(G+ L)

+
sin2 g

(
G2 −K2

)
4G2L(G+ L)

−
(
G2 +K2

)
4G2L2

=
cos(2g)(G−K)(G+K)(G− L)

4G2L2(G+ L)
−
(
G2 +K2

)
4G2L2

=
cos(2g)(G−K)(G+K)(G− L)−

(
G2 +K2

)
(G+ L)

4G2L2(G+ L)
.



24 M. LI, Q. LIU EJDE-2021/63

5.3. Averaged system. Now we state the averaged system of (4.4) into an explicit
expression. Since H10 is independent of g, the averaged system (5.2) becomes

dg

dl
= εL3

(∂H10

∂G
(g,G;L,K) +

∂H11

∂G
(g,G;L,K)

)
:= εF1(g,G;h,K0),

dG

dl
= −εL3 ∂H11

∂g
(g,G;L,K) := εF2(g,G;h,K0),

(5.7)

where K = K0 and L = 1/
√
−2h are considered as constants. Taking the partial

derivative of H11 with respect to the variable g and substituting L = 1/
√
−2h, we

obtain

F2(g,G;h,K0) = −
sin(2g)

(
1−G

√
−2h

) (
G2 −K2

0

)
2
√
−2hG2

(
G
√
−2h+ 1

) .

Moreover, taking the partial derivative of H10 with respect to the variable G and
substituting L = 1/

√
−2h, we obtain

F1(g,G;h,K0) = D(G;h)− 1
√
−2hG3

(
2G
√
−h+

√
2
)2

×
(

cos(2g)
(
−G3

√
−2h+ 2G2hK2

0 +G
√
−2hK2

0 +K2
0

)
+K2

0

(
2G2h− 2G

√
−2h− 1

))
,

where for p > 1,

D(G;h) =L3 ∂H10

∂G
(g,G;L,K)

=− β
√
−2h

2p−3
p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

1

G

×
(
−m

( 1

G
√
−2h

− 1
)m−1

(G
√
−2h)

−p

+ (1− p)
( 1

G
√
−2h

− 1
)m(

G
√
−2h

)1−p)
=− β

√
−2h

2p−3
p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

1

G

×
(

(1− p)(1−G
√
−2h)−m

)( 1

G
√
−2h

− 1
)m−1(

G
√
−2h

)−p
=− β

√
−2h

p−3
G−p−1

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

×
(

(1− p)(1−G
√
−2h)−m

)(1−G
√
−2h

G
√
−2h

)m−1
:=D1(G;h, p),

and for p ≤ 1,

D(G;h) =L3 ∂H10

∂G
(g,G;L,K)



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 25

=− β2p−1
√
−2h

2p−2
1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

(1−G
√
−2h

1 +G
√
−2h

)m−1
×
(
G
√
−2h+ 1

)−p−1(
(−2m− p+ 1)− (1− p)G

√
−2h

)
=− β2p−1

√
−2h

2p−2(
G
√
−2h+ 1

)−p−1 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(−2m− p+ 1)− (1− p)G
√
−2h

)(1−G
√
−2h

1 +G
√
−2h

)m−1
=− β2p−1

√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(1− p)(1−G
√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1
:=D2(G;h, p).

6. Proof of main results

In this section, we give the proofs of Theorems 1.1 and 1.2. The proofs are based
on the method of averaging; For example, see [25, 28, 26, 27]. Generally, averaging
method involves two steps: transforming to standard form; solving the averaging
equation. Now we compute the equilibrium points of the averaged system (5.7).
From the second equation of the averaged system (5.7), it is easy to have that
sin(2g)

(
1−G

√
−2h

) (
G2 −K2

0

)
= 0.

Proof of Theorem 1.1. In this subsection, we always assume that p > 1. We
separate the computation of equilibrium points in three cases.

Case I:
(
1−G

√
−2h

)
= 0. This means that G∗ = 1√

−2h . The equilibrium points

must satisfy the first equation of the averaged system (5.7). Inserting G∗ back into
F1(g,G;h,K0) = 0, we obtain

F1

(
g,

1√
−2h

;h,K0

)
= β2p−2p(p− 1)(−h)p−1 +

1

8
cos(2g)

(
2hK2

0 + 1
)
− hK2

0 = 0.

Therefore,

cos(2g) =
8
(
− β2p−2p(p− 1)(−h)p−1 + hK2

0

)
2hK2

0 + 1

which possesses the solutions

g∗ = ±1

2
arccos

(8
(
−β2p−2p(p− 1)(−h)p−1 + hK2

0

)
2hK2

0 + 1

)
,

which exist when

8
(
−β2p−2p(p− 1)(−h)p−1 + hK2

0

)
2hK2

0 + 1
∈ [−1, 1].

Therefore, we have two solutions

(g∗, G∗;h,K0) =
(
± 1

2
arccos

(8
(
−β2p−2p(p− 1)(−h)p−1 + hK2

0

)
2hK2

0 + 1

)
,

1√
−2h

)
.



26 M. LI, Q. LIU EJDE-2021/63

The Jacobian determinant of the averaged equations at these two solutions is∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=
1

16

(
2hK2

0 + 1
)2 (

1−
4
(
β2pp(p− 1)(−h)p + 4h2K2

0

)2
h2 (2hK2

0 + 1)
2

)
.

Only when

K0 6=
1√
−2h

, |β| 6=
8hK0

2 +
∣∣2hK0

2 + 1
∣∣

2p+1p(p− 1)|h|p−1
,

we can get the nonzero Jacobian determinant.
There exists two 2π-periodic solutions (g(l, ε), G(l, ε)) = (g∗ +O(ε), G∗ +O(ε))

of the system (4.4). In order for k(l) to be 2π-periodic, we plug (g(l, ε), G(l, ε))
back into the second equation of the system (4.2), and obtain

dk

dl
= εL3 ∂H1

∂K
(l, g(l, ε), L,G(l, ε),K) +O

(
ε2
)
.

Next, we prove that there exists K = K0 such that

F3(K0, ε) :=

∫ 2π

0

L3 ∂H1

∂K
(l, g(l, ε), L,G(l, ε),K0)dl +O(ε) = 0.

Since the H1 and ∂H1/∂K are continuous, and we obtain

F3(K0, ε) =

∫ 2π

0

L3 ∂H1

∂K
(l, g(l, ε), L,G(l, ε),K0)dl +O(ε)

=L3 ∂

∂K

∫ 2π

0

H1(l, g(l, ε), L,G(l, ε),K0)dl +O(ε).

Substituting L = 1/
√
−2h and (g(l, ε), G(l, ε)) = (g∗ +O(ε), G∗ +O(ε)), we obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(
1−G∗

√
−2h

)
+G∗

√
−2h+ 1

)
√
−2h(G∗)

2 (
G∗
√
−2h+ 1

) +O(ε)

= −
√
−2hπK0 +O(ε).

Hence there exists (K0, ε) = (0, 0) such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using implicit function theorem we obtain that there exists K0 = K0(ε) such that
F3(K0(ε), ε) = 0. Therefore, there exists K = K0(ε) such that k(l) is 2π-periodic.
We plug K0 = 0 back into the two solutions (g∗, G∗;h,K0), and obtain

(g∗, G∗;h) =
(
± 1

2
arccos

(
−β2p+1p(p− 1)(−h)p−1

)
,

1√
−2h

)
with

|β| < 1

2p+1p(p− 1)|h|p−1
.

Then, statement (i) of Theorem 1.1 is proved.

Case II:
(
G2 −K2

0

)
= 0. This equation have two solutions G∗ = ±K0. Because

the equilibrium points must satisfy the first equation of the averaged system (5.7).
Inserting G∗ back into F1(g,G;h,K0) = 0, we obtain

F1(g,±K0;h,K0) = D1(±K0;h, p) +
1

2
√
−2hK0

∓
cos(2g)

(
1−
√
−2hK0

)
2
√
−2hK0

(
1 +
√
−2hK0

) = 0.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 27

Because we consider the existence of periodic orbits of the system, when K0 near
to origin. K0 = 0 is the singular point of function F1(g,±K0;h,K0), so this case
can not be discussed.

Case III: sin(2g) = 0. This equation has four solutions g∗ = 0,±π2 , π. Next we
consider four subcases.

Subcase A: g∗ = 0. The equilibrium points must satisfy the first equation of the
averaged system (5.7). Inserting g∗ back into F1(g,G;h,K0) = 0, we obtain

F1(0, G;h,K0) =D1(G;h, p)− 1
√
−2hG3

(
2G
√
−h+

√
2
)2

×
((
−G3

√
−2h+ 2G2hK2

0 +G
√
−2hK2

0 +K2
0

)
+K2

0

(
2G2h− 2G

√
−2h− 1

))
= 0.

(6.1)

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain

F1(0, G;h,K0) =D1(G;h, p) +
1

2
(
G
√
−2h+ 1

)2 +O(K0)

=− β
√
−2h

p−3
G−p−1

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

×
(

(1− p)(1−G
√
−2h)−m

)(1−G
√
−2h

G
√
−2h

)m−1
+

1

2
(
G
√
−2h+ 1

)2 +O(K0)

=− β
√
−2h

p−2
G−p

(
(1− p) +

(p− 1)(p− 2)

2

×
(
− p1−G

√
−2h

G
√
−2h

− 1
) p−2∑
m=2

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

×
(

(1− p)(1−G
√
−2h)−m

)(1−G
√
−2h

G
√
−2h

)m−1)
+

1

2
(
G
√
−2h+ 1

)2 +O(K0).

(6.2)

From this equality, we obtain that F1(0, G;h,K0) have two singularities G =
− 1√
−2h , 0. Therefore G ∈

(
− 1/

√
−2h, 0

)
∪
(
0, 1/
√
−2h

)
.

When G → −1/
√
−2h , we obtain F1(0, G;h,K0) → +∞. When G → 0 and

β > 0, we obtain F1(0, G;h,K0)→ +∞, while if β < 0 we obtain F1(0, G;h,K0)→
−∞. When G→ 1/

√
−2h and

β > − 1

2p+1(−h)p−1p(p− 1)
,

we have F1(0, G;h,K0) > 0, while if

β < − 1

2p+1(−h)p−1p(p− 1)

we have F1(0, G;h,K0) < 0.



28 M. LI, Q. LIU EJDE-2021/63

Using zero point theorem, we obtain that when

β ∈
(
− 1

2p+1(−h)p−1p(p− 1)
, 0
)
,

there exist two solutions G̃1(K0, h) ∈ (−1/
√
−2h, 0) and G̃2(K0, h) ∈ (0, 1/

√
−2h)

of (6.1). Therefore, when

β ∈
(
− 1

2p+1(−h)p−1p(p− 1)
, 0
)
,

we have two solutions (g∗, G∗(1);K0, h) = (0, G̃1(K0, h)) and (g∗, G∗(2);K0, h) =

(0, G̃2(K0, h)).
Next, we calculate the Jacobian determinant of the averaged equation at these

solutions. Firstly, we consider the derivative of D1(G;h, p) with respect to G,

D′1(G;h, p) =
d

dG
D1(G;h, p)

=− β
√
−2h

p−4
G−p−3

p−2∑
m=0

2−mCmp−2
(m+ p− 2)!

m!(p− 2)!

(1−G
√
−2h

G
√
−2h

)m−2
×
((

(1− p)
(
1−G

√
−2h

)
−m

)(
(−1− p)

(
1−G

√
−2h

)
+ 1−m

)
+G(p− 1)

√
−2h

(
1−G

√
−2h

))
.

Therefore we can obtain the derivative of F1(g,G;h,K0) with respect to G. Then

we plug (g∗, G∗(i);K0, h) = (0, G̃i(K0, h)), i = 1, 2 back into the derivative, and
obtain

d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

= D′1(G̃i;h, p)−
(√
−2hG̃3

i − 6G̃2
ihK0

2 + 3
√
−2hG̃iK0

2 +K0
2
)

G̃3
i

(√
−2hG̃i + 1

)3
=
(
D′1(G̃i;h, p)G̃

3
i

(√
−2hG̃i + 1

)3
−
(√
−2hG̃3

i − 6G̃2
ihK0

2 + 3
√
−2hG̃iK0

2 +K0
2
))/(

G̃3
i

(√
−2hG̃i + 1

)3)
.

The Jacobian determinant of the averaged equations with these solutions are

∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗(i))

= − d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

· d

dg
F2(g,G;h,K0)

∣∣
(g∗,G∗(i))

=

(√
−2hG̃i − 1

)(
G̃2
i −K0

2
)

√
−2hG̃2

i

(
1 +
√
−2hG̃i

) ( d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

)
.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 29

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗(i))

=

(√
−2hG̃i − 1

)
G̃2
i√

−2hG̃2
i

(
1 +
√
−2hG̃i

)( d

dG
F1(g,G;h)

∣∣
(g∗,G∗(i))

)
+O(K0)

=

(√
−2hG̃i − 1

)
√
−2h

(
1 +
√
−2hG̃i

) (D′1(G̃i;h, p)
(√
−2hG̃i + 1

)3 −√−2h
)(√

−2hG̃i + 1
)3 +O(K0)

=

(√
−2hG̃i − 1

)(
D′1(G̃i;h, p)

(√
−2hG̃i + 1

)3 −√−2h
)

√
−2h

(√
−2hG̃i + 1

)4 +O(K0).

(6.3)

Because G̃i 6= ±1/
√
−2h, to get a nonzero Jacobian determinant, expression (6.3)

needs to satisfy (
D′1(G̃i;h, p)

(√
−2hG̃i + 1

)3 −√−2h
)
6= 0.

Because G̃i is the nondegenerate zero of (6.2) and

d

dG

(
D1(G;h, p) +

1

2
(
G
√
−2h+ 1

)2 ) = D′1(G;h, p)−
√
−2h(

G
√
−2h.+ 1

)3 ,
we have(

D′1(G̃i;h, p)
(√
−2hG̃i + 1

)3 −√−2h
)

= D′1(G̃i;h, p)−
√
−2h(

G̃i
√
−2h+ 1

)3 6= 0.

Then there exists 2π-periodic solutions (g(l, ε), G(l, ε)) =
(
g∗ +O(ε), G∗(i) +O(ε)

)
,

i = 1, 2 of the system (4.4). Similarly, in order for k(l) to be 2π-periodic, we plug
(g(l, ε), G(l, ε)) back into F3(K0, ε), and obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗(i)

√
−h
)

+ 2G∗(i)
√
−h+

√
2
)

√
−2h(G∗(i))

2 (
2G∗(i)

√
−h+

√
2
) +O(ε)

= − 2πK0

G̃i
(
G̃i
√
−2h+ 1

) +O(ε).

There exists (K0, ε) = (0, 0), such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using implicit function theorem we obtain that there exists K0 = K0(ε), such that
F3(K0(ε), ε) = 0. Therefore there exists K = K0(ε) such that k(l) is 2π-periodic.

Then, statement (ii) of Theorem 1.1 is proved when g∗ = 0.

Subcase B: g∗ = π. The proof of statement (ii) of Theorem 1.1 when g∗ = π is
completely similar to case g∗ = 0.



30 M. LI, Q. LIU EJDE-2021/63

Subcase C: g∗ = π
2 . The equilibrium points must satisfy the first equation of the

averaged system (5.7). Inserting g∗ back into F1(g,G;h,K0) = 0, we obtain

F1

(π
2
, G;h,K0

)
=D1(G;h, p) +

1
√
−2hG3

(
2G
√
−h+

√
2
)2

×
((
−G3

√
−2h+ 2G2hK2

0 +G
√
−2hK2

0 +K2
0

)
−K2

0

(
2G2h− 2G

√
−2h− 1

))
= 0.

(6.4)

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain

F1

(π
2
, G;h,K0

)
= D1(G;h, p)− 1

2
(
G
√
−2h+ 1

)2 +O(K0)

= −β
√
−2h

p−3
G−p−1

p−2∑
m=0

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

×
(

(1− p)(1−G
√
−2h)−m

)(1−G
√
−2h

G
√
−2h

)m−1
− 1

2
(
G
√
−2h+ 1

)2 +O(K0)

= −β
√
−2h

p−2
G−p

(
(1− p) +

(p− 1)(p− 2)

2

(
− p1−G

√
−2h

G
√
−2h

− 1
)

+

p−2∑
m=2

Cmp−22−m
(m+ p− 2)!

m!(p− 2)!

(
(1− p)(1−G

√
−2h)−m

)
×
(1−G

√
−2h

G
√
−2h

)m−1)
− 1

2
(
G
√
−2h+ 1

)2 +O(K0).

(6.5)

From this equality, we obtain that F1(π2 , G;h,K0) have two singularities G =

−1/
√
−2h, 0. Therefore G ∈ (−1/

√
−2h, 0) ∪ (0, 1/

√
−2h).

When G → −1/
√
−2h, we obtain F1(π/2, G;h,K0) → −∞. When G → 0 and

β > 0, we obtain F1(π/2, G;h,K0) → +∞. When G → 0 and β < 0, we obtain
F1(π/2, G;h,K0)→ −∞. When G→ 1/

√
−2h, and

β >
1

2p+1(−h)p−1p(p− 1)
,

we obtain F1(π/2, G;h,K0) > 0.
Similarity, when G→ 1/

√
−2h, and

β <
1

2p+1(−h)p−1p(p− 1)

we obtain F1(π/2, G;h,K0) < 0.
Using the zero point theorem, we obtain that when

β ∈
(

0,
1

2p+1(−h)p−1p(p− 1)

)
,



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 31

there exist two solutions G1(K0, h) ∈ (−1/
√
−2h, 0) and G2(K0, h) ∈ (0, 1/

√
−2h)

of (6.4). Therefore, when

β ∈
(

0,
1

2p+1(−h)p−1p(p− 1)

)
,

we have two solutions

(g∗, G∗(1);K0, h) =
(π

2
, G1(K0, h)

)
,

(g∗, G∗(2);K0, h) =
(π

2
, G2(K0, h)

)
.

Next, we calculate the Jacobian determinant of the averaged equation with these
solutions. Firstly, we consider the derivative of F1(g,G;h,K0) with respect to G.
Then we insert

(g∗, G∗(i);h,K0) =
(π

2
, Gi(h,K0)

)
, i = 1, 2

back into the derivative, and obtain

d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

= D′1(Gi;h, p) +
−2hG

4

i + 12hG
2

iK0
2 −Gi8

√
−2hK0

2 − 3K0
2

G
4

i

√
−2h

(
Gi
√
−2h+ 1

)3
=

(
D′1(Gi;h, p)G

4

i

√
−2h

(
Gi
√
−2h+ 1

)3
+
(
− 2hG

4

i + 12hG
2

iK0
2 −Gi8

√
−2hK0

2 − 3K0
2
))

÷
(
G

4

i

√
−2h

(
Gi
√
−2h+ 1

)3 )
.

The Jacobian determinant of the averaged equations with these solutions are

∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗(i))

= − d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

· d
dg
F2(g,G;h,K0)

∣∣
(g∗,G∗(i))

=

(
Gi
√
−2h− 1

) (
G

2

i −K0
2
)

−G2

i

√
−2h

(
1 +Gi

√
−2h

) ( d

dG
F1(g,G;h,K0)

∣∣
(g∗,G∗(i))

)



32 M. LI, Q. LIU EJDE-2021/63

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗(i))

=

(
Gi
√
−2h− 1

)
G

2

i

−G2

i

√
−2h

(
1 +
√
−2hGi

)( d

dG
F1(g,G;h)

∣∣
(g∗,G∗(i))

)
+O(K0)

=

(
Gi
√
−2h− 1

)
−
√
−2h

(
1 +Gi

√
−2h

) (D′1(Gi;h, p)
(
Gi
√
−2h+ 1

)3
+
√
−2h

)(
Gi
√
−2h+ 1

)3
+O(K0)

=

(
Gi
√
−2h− 1

)(
D′1(Gi;h, p)

(
Gi
√
−2h+ 1

)3
+
√
−2h

)
−
√
−2h

(
Gi
√
−2h+ 1

)4 +O(K0).

(6.6)

Because Gi 6= ±1/
√
−2h, to obtain the nonzero Jacobian determinant, (6.6) needs

to satisfy (
D′1(Gi;h)

(
Gi
√
−2h+ 1

)3
+
√
−2h

)
6= 0.

Because Gi is the nondegenerate zero of (6.5) and

d

dG

(
D1(G;h, p)− 1

2
(
G
√
−2h+ 1

)2) = D′1(G;h, p) +

√
−2h(

G
√
−2h+ 1

)3 ,
we have(

D′1(Gi;h, p)
(√
−2hGi + 1

)3
+
√
−2h

)
= D′1(Gi;h, p) +

√
−2h(

Gi
√
−2h+ 1

)3 6= 0.

There exists 2π-periodic solutions (g(l, ε), G(l, ε)) =
(
g∗ +O(ε), G∗(i) +O(ε)

)
, i =

1, 2 of the system (4.4). Similarly, in order for k(l) to be 2π-periodic, we plug
(g(l, ε), G(l, ε)) back into F3(K0, ε) to obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗(i)

√
−h
)

+ 2G∗(i)
√
−h+

√
2
)

√
−2h(G∗(i))

2 (
2G∗(i)

√
−h+

√
2
) +O(ε)

=
2πK0

Gi
(
Gi
√
−2h+ 1

) +O(ε).

There exists (K0, ε) = (0, 0) such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using the implicit function theorem we obtain that there exists K0 = K0(ε), such
that F3(K0(ε), ε) = 0. Therefore there exists K = K0(ε) such that k(l) is 2π-
periodic. Then, statement (iii) of Theorem 1.1 is proved when g∗ = π/2.

Subcase D: g∗ = −π2 . The proof of statement (iii) of Theorem 1.1 with g∗ = −π/2
is completely similar to that of case g∗ = π/2.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 33

6.1. Proof of Theorem 1.2. In this subsection, we always assume that p ≤ 1.
We separate the computation of the equilibrium points in three cases.

Case I:
(
1−G

√
−2h

)
= 0. This means that G∗ = 1/

√
−2h. The equilibrium

points must satisfy the first equation of the averaged system (5.7). Inserting G∗

back into F1(g,G;h,K0) = 0, we obtain

F1

(
g,

1√
−2h

;h,K0

)
= −β2p−2(1−p)p(−h)p−1+

1

8
cos(2g)

(
2hK0

2 + 1
)
−hK0

2 = 0.

Therefore,

cos(2g) =
8
(
hK0

2 + β2p−2(1− p)p(−h)p−1
)

2hK0
2 + 1

.

It possesses the solutions

g∗ = ±1

2
arccos

(8
(
hK0

2 + β2p−2(1− p)p(−h)p−1
)

2hK0
2 + 1

)
,

which exist when

8
(
hK0

2 + β2p−2(1− p)p(−h)p−1
)

2hK0
2 + 1

∈ [−1, 1].

Therefore, we have two solutions:

(g∗, G∗;h,K0) =
(
± 1

2
arccos

(8
(
hK2 + β2p−2(1− p)p(−h)p−1

)
2hK2 + 1

)
,

1√
−2h

)
.

The Jacobian determinant of the averaged equations at these two solutions is∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=
1

16

(
2hK2

0 + 1
)2 (

1−
4
(
4h2K2

0 + β2p(p− 1)p(−h)p
)2

h2 (2hK2
0 + 1)

2

)
.

Only when

K0 6=
1√
−2h

, |β| 6=
8hK0

2 +
∣∣2hK0

2 + 1
∣∣

2p+1(p− 1)p|h|p−1
we can get the nonzero Jacobian determinant.

There exist two 2π-periodic solutions (g(l, ε), G(l, ε)) = (g∗ +O(ε), G∗ +O(ε)) of
system (4.4). Similarly, in order for k(l) to be 2π-periodic, we plug (g(l, ε), G(l, ε))
back into F3(K, ε), and obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗

√
−h
)

+ 2G∗
√
−h+

√
2
)

√
−2h(G∗)

2 (
2G∗
√
−h+

√
2
) +O(ε)

= −
√
−2hπK0 +O(ε).

There exists (K0, ε) = (0, 0) such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using implicit function theorem we obtain that there exists K0 = K0(ε) such that
F3(K0(ε), ε) = 0. Therefore there exists K0 = K0(ε) such that k(l) is 2π-periodic.
We plug K0 = 0 back into the two solutions (g∗, G∗;h,K0) and obtain

(g∗, G∗;h) =
(
± 1

2
arccos

(
β2p+1(1− p)p(−h)p−1

)
,

1√
−2h

)
,



34 M. LI, Q. LIU EJDE-2021/63

with

|β| < 1

2p+1(p− 1)p|h|p−1
.

Then, statement (i) of Theorem 1.2 is proved.

Case II:
(
G2 −K2

0

)
= 0. This equation have two solutions G∗ = ±K0. Due to

the equilibrium points must satisfy the first equation of the averaged system (5.7).
Inserting G∗ back into F1(g,G;h,K0) = 0, we obtain

F1(g,±K0;h,K0) = D2(±K0;h, p) +
1

2
√
−2hK0

∓
cos(2g)

(
1−
√
−2hK0

)
2
√
−2hK0

(
1 +
√
−2hK0

) = 0.

We consider the existence of periodic orbits of the system for the case that K0 is
near the origin. Since K0 = 0 is the singular point of function F1(g,±K0;h,K0),
this case is not feasible.

Case III: sin(2g) = 0. We have four solutions g∗ = 0,±π2 , π. Next consider four
subcases.

Subcase A: g∗ = 0. The equilibrium points must satisfy the first equation of the
averaged system (5.7). Inserting g∗ back into F1(g,G;h,K0) = 0, we obtain

F1(0, G;h,K0) =D2(G;h, p)− 1
√
−2hG3

(
2G
√
−h+

√
2
)2

×
((
−G3

√
−2h+ 2G2hK2

0 +G
√
−2hK2

0 +K2
0

)
+K2

0

(
2G2h− 2G

√
−2h− 1

))
= 0.

(6.7)

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain

F1(0, G;h,K0)

= D2(G;h, p) +
1

2
(
G
√
−2h+ 1

)2 +O(K0)

= −β2p−1
√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(1− p)(1−G
√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1
+

1

2
(
G
√
−2h+ 1

)2 +O(K0)

= −β2p−1
√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 (
(1− p)(1 +G

√
−2h)

+ (1− p)2
(

(1− p)(1−G
√
−2h)− 2

)
+

1−p∑
m=2

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(1− p)(1−G
√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1)
+

1

2
(
G
√
−2h+ 1

)2 +O(K0).

(6.8)



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 35

When p = 1 or p = 0, we obtain

F1(0, G;h,K0) =
1

2
(
G
√
−2h+ 1

)2 +O(K0),

there is no solutionG(K0, h). When p < 0, from (6.8), we obtain that F1(0, G;h,K0)
has one singularity G = − 1√

−2h . Therefore G ∈
(
− 1/

√
−2h, 1/

√
−2h

)
.

When G→ 0, we obtain F1(0, G;h,K0)→ 1/2. When G→ 1/
√
−2h and

β > − 1

2p+1(p− 1)p(−h)p−1
,

we obtain F1(0, G;h,K0) > 0, while if

β < − 1

2p+1(p− 1)p(−h)p−1
,

we have F1(0, G;h,K0) < 0.
Using the zero point theorem, we obtain that when

β ∈
(
−∞,− 1

2p+1(p− 1)p(−h)p−1

)
,

there exists a solution G̃(K0, h) ∈
(

0, 1√
−2h

)
of (6.7).

Therefore, when p < 0 and

β ∈
(
−∞,− 1

2p+1(p− 1)p(−h)p−1

)
we have a solution (g∗, G∗;K0, h) = (0, G̃(K0, h)).

Next, we calculate the Jacobian determinant of the averaged equation at this
solution. Firstly, we consider the derivative of D2(G;h, p) with respect to G,

D′2(G;h, p)

=
d

dG
D2(G;h, p)

= −β2p−1
√
−2h

2p−2(
G
√
−2h+ 1

)−p−3 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(1−G

√
−2h

1 +G
√
−2h

)m−2((√
−2h(−p− 1)

(
1−G

√
−2h

)
+ (1−m)(1 +

√
−2h)

)
×
(

(1− p)
(
1−G

√
−2h

)
− 2m

)
−G
√
−2h(1− p)

(
1−G

√
−2h

))
.

Therefore we can obtain the derivative of F1(g,G;h,K0) with respect to G. Then

we plug (g∗, G∗;K0, h) = (0, G̃(K0, h)) back into the derivative, to obtain

d

dG
F1(g,G;h,K0)|(g∗,G∗)

= D′2(G̃;h, p)−
(√
−2hG̃3 − 6G̃2hK0

2 + 3
√
−2hG̃K0

2 +K0
2
)

G̃3
(√
−2hG̃+ 1

)3
=
D′2(G̃;h, p)G̃3

(√
−2hG̃+ 1

)3
−
(√
−2hG̃3 − 6G̃2hK0

2 + 3
√
−2hG̃K0

2 +K0
2
)

G̃3
(√
−2hG̃+ 1

)3 .



36 M. LI, Q. LIU EJDE-2021/63

The Jacobian determinant of the averaged equations with this solution is∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=− d

dG
F1(g,G;h,K0)|(g∗,G∗) ·

d

dG
F1(g,G;h,K0)|(g∗,G∗)

=

(√
−2hG̃− 1

)(
G̃2 −K0

2
)

√
−2hG̃2

(
1 +
√
−2hG̃

) ( d

dG
F1(g,G;h,K0)|(g∗,G∗)

)
.

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=

(√
−2hG̃− 1

)
G̃2

√
−2hG̃2

(
1 +
√
−2hG̃

)( d

dG
F1(g,G;h)|(g∗,G∗)

)
+O(K0)

=

(√
−2hG̃− 1

)
√
−2h

(
1 +
√
−2hG̃

)
(
D′2(G̃;h, p)

(√
−2hG̃+ 1

)3 −√−2h
)

(√
−2hG̃+ 1

)3 +O(K0)

=

(√
−2hG̃− 1

)(
D′2(G̃;h, p)

(√
−2hG̃+ 1

)3 −√−2h
)

√
−2h

(√
−2hG̃+ 1

)4 +O(K0).

(6.9)

Because G̃ 6= ±1/
√
−2h, to obtain the nonzero Jacobian determinant, (6.9) needs

to satisfy (
D′2(G̃;h, p)

(√
−2hG̃+ 1

)3 −√−2h
)
6= 0.

Because G̃ is the nondegenerate zero of (6.8) and

d

dG

(
D2(G;h, p) +

1

2
(
G
√
−2h+ 1

)2) = D′2(G;h, p)−
√
−2h(

G
√
−2h.+ 1

)3 .
Therefore,(

D′2(G̃;h, p)
(√
−2hG̃+ 1

)3 −√−2h
)

= D′2(G̃;h, p)−
√
−2h(

G̃
√
−2h+ 1

)3 6= 0.

There exists 2π-periodic solutions (g(l, ε), G(l, ε)) = (g∗ +O(ε), G∗ +O(ε)) of the
system (4.4). Similarly, in order for k(l) to be 2π-periodic, we plug (g(l, ε), G(l, ε))
back into F3(K0, ε), and obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗

√
−h
)

+ 2G∗
√
−h+

√
2
)

√
−2h(G∗)

2 (
2G∗
√
−h+

√
2
) +O(ε)

= − 2πK0

G̃
(
G̃
√
−2h+ 1

) +O(ε).

There exists (K0, ε) = (0, 0) such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using implicit function theorem we obtain that there exists K0 = K0(ε), such that
F3(K0(ε), ε) = 0. Therefore there exists K = K0(ε) such that k(l) is 2π-periodic.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 37

When p = 1 or p = 0, we obtain D2(G;h, p) = 0. Solving F1(0, G;h,K0) = 0
with respect to G, we obtain

G∗± = −
√
−2hK2

0 ±
√
−K2

0 − 2hK4
0 .

Then, we plug g∗ and G∗± back into F3(K0, ε) to obtain

F±3 (K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗+

√
−h
)

+ 2G∗+
√
−h+

√
2
)

√
−2h(G∗+)

2 (
2G∗+
√
−h+

√
2
) +O(ε)

= −
2
√
−2hπK0

(√
hK2

0 (1 + 2hK2
0 )∓

√
2hK2

0

)
√

2
√
hK2

0 (1 + 2hK2
0 )
(√

hK2
0 (1 + 2hK2

0 )±
√

2hK2
0

) +O(ε).

There not exists K0(ε) such that F3(K0(ε), ε) = 0. Therefore there not exists
K = K0(ε) such that k(l) is 2π-periodic. Then, statement (ii) of Theorem 1.2 is
proved.

Subcase B: g∗ = π. The proof of statement (ii) of Theorem 1.2 when g∗ = π is
completely similar to case g∗ = 0.

Subcase C: g∗ = π
2 . The equilibrium points must satisfy the first equation of the

averaged system (5.7). Inserting g∗ back into F1(g,G;h,K0) = 0, we obtain

F1

(π
2
, G;h,K0

)
=D2(G;h, p) +

1
√
−2hG3

(
2G
√
−h+

√
2
)2

×
((
−G3

√
−2h+ 2G2h K2

0 +G
√
−2hK2

0 −K2
0

)
−K2

0

(
2G2h− 2G

√
−2h− 1

))
= 0.

(6.10)

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain

F1

(π
2
, G;h,K0

)
=D2(G;h, p)− 1(

2G
√
−h+

√
2
)2 +O(K0)

=− β2p−1
√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 1−p∑
m=0

Cm1−p
(1− p)!

m!(−m− p+ 1)!

×
(

(1− p)(1−G
√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1
− 1

2
(
G
√
−2h+ 1

)2 +O(K0)

=− β2p−1
√
−2h

2p−2 (
G
√
−2h+ 1

)−p−1 (
(1− p)(1 +G

√
−2h)

+ (1− p)2
(

(1− p)(1−G
√
−2h)− 2

)
+

1−p∑
m=2

Cm1−p
(1− p)!

m!(−m− p+ 1)!(
(1− p)(1−G

√
−2h)− 2m

)(1−G
√
−2h

1 +G
√
−2h

)m−1)



38 M. LI, Q. LIU EJDE-2021/63

− 1

2
(
G
√
−2h+ 1

)2 +O(K0). (6.11)

When p = 1 or p = 0, we obtain

F1

(π
2
, G;h,K0

)
= − 1

2
(
G
√
−2h+ 1

)2 +O(K0),

there is no solutionG(K0, h). When p < 0, from (6.11), we obtain that F1(π2 , G;h,K0)

has one singularity G = − 1√
−2h . Therefore G ∈

(
− 1/

√
−2h, 1/

√
−2h

)
.

When G→ 0, we obtain F1(π/2, G;h,K0)→ − 1
2 . When G→ 1√

−2h and

β >
1

2p+1(p− 1)p(−h)p−1
,

we obtain F1 (π/2, G;h,K0) > 0, while if

β <
1

2p+1(p− 1)p(−h)p−1
,

we have F1 (π/2, G;h,K0) < 0.
Using the zero point theorem, we obtain that when

β ∈
( 1

2p+1(p− 1)p(−h)p−1
,+∞

)
,

there exists a solution G(K0, h) ∈ (0, 1/
√
−2h) of (6.10). Therefore, when p < 0

and

β ∈
(

0,
1

2p+1(p− 1)p(−h)p−1

)
,

we have a solution

(g∗, G∗;K0, h) =
(π

2
, G(K0, h)

)
.

Next, we calculate the Jacobian determinant of the averaged equation at this
solution. Firstly, we consider the derivative of F1(g,G;h,K0) with respect to G.
Then we put (g∗, G∗;h,K0) = (π2 , G(h,K0)) back into the derivative, and obtain

d

dG
F1(g,G;h,K0)|(g∗,G∗)

= D′2(G;h, p) +
−2hG

4
+ 12hG

2
K0

2 − 8G
√
−2hK0

2 − 3K0
2

G4
√
−2h

(
G
√
−2h+ 1

)3
=

(
D′2(G;h, p)G

4√−2h
(
G
√
−2h+ 1

)3
+
(
−2hG

4
+ 12hG

2
K0

2 − 8G
√
−2hK0

2 − 3K0
2
))

÷
(
G

4√−2h
(
G
√
−2h+ 1

)3)
The Jacobian determinant of the averaged equation with this solution is∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=− d

dG
F1(g,G;h,K0)|(g∗,G∗) ·

d

dg
F2(g,G;h,K0)|(g∗,G∗)

=

(
G
√
−2h− 1

) (
G

2 −K0
2
)

−G2√−2h
(
1 +G

√
−2h

) ( d

dG
F1(g,G;h,K0)|(g∗,G∗)

)
.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 39

We restrict K0 to a sufficiently small neighborhood of the origin, and obtain∣∣∣∂(F1, F2)

∂(g,G)

∣∣∣
(g∗,G∗)

=

(
G
√
−2h− 1

)
G

2

−G2√−2h
(
1 +G

√
−2h

)( d

dG
F1(g,G;h)|(g∗,G∗)

)
+O(K0)

=

(
G
√
−2h− 1

)
−
√
−2h

(
1 +G

√
−2h

) (D′2(G;h, p)
(
G
√
−2h+ 1

)3
+
√
−2h

)(
G
√
−2h+ 1

)3
+O(K0)

=

(
G
√
−2h− 1

) (
D′2(G;h, p)

(
G
√
−2h+ 1

)3
+
√
−2h

)
−
√
−2h

(
G
√
−2h+ 1

)4 +O(K0).

(6.12)

Because G 6= ±1/
√
−2h, to obtain the nonzero Jacobian determinant, (6.12) needs

to satisfy (
D′2(G;h, p)

(
G
√
−2h+ 1

)3
+
√
−2h

)
6= 0.

Because G is the nondegenerate zero of (6.11) and

d

dG

(
D2(G;h, p)− 1

2
(
G
√
−2h+ 1

)2) = D′2(G;h, p) +

√
−2h(

G
√
−2h+ 1

)3 .
Therefore,(

D′2(G;h, p)
(√
−2hG+ 1

)3
+
√
−2h

)
= D′2(G;h, p) +

√
−2h(

G
√
−2h+ 1

)3 6= 0.

There exists a 2π-periodic solution (g(l, ε), G(l, ε)) = (g∗ +O(ε), G∗ +O(ε)) of sys-
tem (4.4). Similarly, in order for k(l) to be 2π-periodic, we plug (g(l, ε), G(l, ε))
back into F3(K0, ε), obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗

√
−h
)

+ 2G∗
√
−h+

√
2
)

√
−2h(G∗)

2 (
2G∗
√
−h+

√
2
) +O(ε)

=
2πK0

G
(
G
√
−2h+ 1

) +O(ε).

There exists (K0, ε) = (0, 0) such that

F3(0, 0) = 0,
d

dK0
F3(K0, ε)|(0,0) 6= 0.

Using the implicit function theorem we obtain that there exists K0 = K0(ε), such
that F3(K0(ε), ε) = 0. Therefore there exists K = K0(ε) such that k(l) is 2π-
periodic.

When p = 1 or p = 0, we obtain D2(G;h, p) = 0. Solving F1(π/2, G;h,K0) = 0
with respect to G, we obtain

G∗ =

(
−h
(√

2hK6 +K4 +K2
))2/3 − 21/6hK2

21/6(−h)5/6
(√

2hK6 +K4 +K2
)1/3 .



40 M. LI, Q. LIU EJDE-2021/63

Then, we plug g∗ and G∗ back into F3(K0, ε), and obtain

F3(K0, ε) = −
πK0

(
− cos(2g∗)

(√
2− 2G∗

√
−h
)

+ 2G∗
√
−h+

√
2
)

√
−2h(G∗)

2 (
2G∗
√
−h+

√
2
) +O(ε)

=
25/6
√
−hhK

(√
2hK6 +K4 +K2

)( (
−h
(√

2hK6 +K4 +K2
))2/3 − 21/6hK2

)2
×
(
−
√

2hK2 +
(
−h
(√

2hK6 +K4 +K2
))1/3

+ 21/3
(
−h
(√

2hK6 +K4 +K2
))2/3)−1

.

There is no K0(ε) such that F3(K0(ε), ε) = 0. Therefore there is not K = K0(ε)
such that k(l) is 2π-periodic. Then, statement (iii) of Theorem 1.2 is proved when
g∗ = π/2.

Subcase D: g∗ = −π/2. The proof of statement (iii) of Theorem 1.2 with g∗ =
−π/2 is completely similar to case g∗ = π/2, and we do not repeat again. Therefore,
we conclude the proof of Theorem 1.2.

7. Conclusion

In this article, the periodic orbits of the spatial anisotropic Kepler problem with
anisotropic perturbations on the energy surface are discussed by using canonical
transformations, Residue theory in complex analysis and the averaging method, see
Theorem 1.1 and Theorem 1.2. Moreover, we also give some applications related
to spatial anisotropic problem. Comparing the previous work, we are not limited
to the study of the spatial anisotropic Manev problem when p = 2, but we extend
p to the integer domain. This is a new result in this context. From the view point
of a physical application, it might be reasonable to use the averaging principle to
replace a mathematical model by the corresponding averaged system, to use the
averaged system to make a prediction. The study of this paper exactly gives a frame
for application of the averaging method to spatial anisotropic Kepler problem with
anisotropic perturbations.

Acknowledgments. This research was supported by the National Natural Science
Foundation (No. 11771105), by the Guangxi Natural Science Foundation (Nos.
2017GXNSFFA198012, 2018GXNSFAA138177), and by the Innovation Project of
Guang Xi Graduate Education (No. YCSW2020157).

References

[1] E. I. Abouelmagd, J. Llibre, L. G. Guirao; Periodic orbits of the planar anisotropic Kepler

problem, International Journal of Bifurcation and Chaos 27 (2017), no. 03, 1–6.

[2] A. Alberti, C. Vidal; New families of symmetric periodic solutions of the spatial anisotropic
Manev problem, Journal of Mathematical Physics 56 (2015), no. 012901, 1–18.

[3] J. Casasayas , J. Llibre; Qualitative analysis of the anisotropic Kepler problem, Memoirs of
the American Mathematical Society 52 (1984), 115.

[4] G. Contopoulos and M. Harsoula; Stability and instability in the anisotropic Kepler problem,

Journal of Physics A General Physics 38 (2005), no. 41, 8897–8920.
[5] B. Cordani; The kepler problem in: geography of order and chaos in mechanics, Birkhauser,

New York, 2003.



EJDE-2021/63 SPATIAL ANISOTROPIC KEPLER PROBLEM 41

[6] S. Craig, F. Diacu, E. A. Lacomba, E. Perez; On the anisotropic Manev problem, Journal of

Mathematical Physics 40 (1999), no. 3, 1359–1375.

[7] J. Delgado, F. N. Diacu, E. A. Lacomba, A. Mingarelli, V. Mioc, E. Perez, C. Stoica; The
global flow of the Manev problem, Journal of Mathematical Physics 37 (1996), no. 6, 2748–

2761.

[8] R. L. Devaney; Collision orbits in the anisotropic Kepler problem, Inventiones Mathematicae
45 (1978), no. 3, 221–251.

[9] R. L. Devaney; Singularities in classical mechanical systems, American Mathematical

Monthly 10 (1981), no. 89, 211–333.
[10] R. L. Devaney; Blowing up singularities in classical mechanical systems, American Mathe-

matical Monthly 89 (1982), no. 8, 535–552.

[11] F. Diacu; Stability in the anisotropic Manev problem, Journal of Physics A General Physics
33 (2000), no. 37, 6573–6578.

[12] F. Diacu, P. Holmes; Celestial encounters: the origins of chaos and stability, Princeton
University Press, 1996.

[13] F. Diacu, E. Pérez-Chavela, M. Santoprete; The kepler problem with anisotropic perturba-

tions, Journal of Mathematical Physics 46 (2005), no. 72701, 1–21.
[14] F. Diacu, E. Pérez-Chavela, M. Santoprete; Saari’s conjecture for the collinear n-body prob-

lem, Transactions of the American Mathematical Society 357 (2009), no. 10, 4215–4223.

[15] F. Diacu, M. Santoprete; Nonintegrability and chaos in the anisotropic Manev problem,
Physica D Nonlinear Phenomena 156 (2001), no. 1, 39–52.

[16] F. Diacu, M. Santoprete; On the global dynamics of the anisotropic Manev problem, Physica

D-nonlinear Phenomena 194 (2004), no. 1-2, 75–94.
[17] A. Einstein; On the quantum theorem of sommerfeld and epstein, Verh Dtsch Ges Pathol 19

(1917), 82–92.

[18] A. H. Escalona-Buend́ıa, E. Pérez-Chavela; Periodic orbits for anisotropic perturbations of
the Kepler problem, Nonlinear Analysis 68 (2008), 591–601.

[19] L. G. Guirao, J. Llibre, J. A. Vera; Periodic orbits of Hamiltonian systems: applications to
perturbed Kepler problems, Chaos Solitons and Fractals 57 (2013), 105–111.

[20] M. C. Gutzwiller; Periodic orbits and classical quantization conditions, Journal of Mathe-

matical Physics 12 (1971), 343–358.
[21] ; The anisotropic Kepler problem in two dimensions, Journal of Mathematical Physics

14 (1973), no. 1, 139–152.

[22] M. C. Gutzwiller; Bernoulli sequences and trajectories in the anisotropic Kepler problem,
Journal of Mathematical Physics 14 (1977), 806–823.

[23] M. C. Gutzwiller; Periodic orbits in the anisotropic Kepler problem, Classical Mechanics and

Dynamical Systems 70 (1981), 69–90.
[24] M. C. Gutzwiller; The quantization of a classically ergodic system, Physica D-Nonlinear

Phenomena 5 (1982), no. 2-3, 183–207.
[25] Z. Li, Q. Liu, K. Zhang; Periodic solutions of nonlinear differential systems by the method

of averaging, Applications of Mathematics 65 (2020), no. 4, 511–542.

[26] Q. Liu, L. Cai; Averaging methods for nonlinear systems with a small parameter via reduction
and topological degree, Nonlinear Analysis: Real World Applications 45 (2019), 461–471.

[27] Q. Liu, L. Huang, G. Jiang; Periodic oscillations of the relativistic pendulum with friction,

Electron. J. Differential Equations 2017 (2017), no. 40, 1–10.
[28] Q. Liu, W. Liu, P. J. Torres, W. Huang; Periodic dynamics of a derivative nonlinear

schrödinger equation with variable coefficients, Applicable Analysis 99 (2020), no. 3, 407–427.
[29] J. Llibre, A. Makhlouf; Periodic orbits of the spatial anisotropic Manev problem, Journal of

Mathematical Physics 53 (2012), no. 12, 312–8920.

[30] J. Llibre, A. E. Teruel, C. Valls, D. L. Fuente; Phase portraits of the two-body problem with

Manev potential, Journal of Physics A General Physics 34 (2001), no. 9, 1919–1934.
[31] J. Llibre, C. Valls; Periodic orbits of the planar anisotropic generalized Kepler problem,

Journal of Mathematical Physics 60 (2019), no. 042901, 1–5.
[32] J. Llibre,P. Yuan; Periodic orbits of the planar anisotropic Manev problem and of the per-

turbed hydrogen atom problem, Qualitative Theory of Dynamical Systems 18 (2019), 969–986.

[33] Miguel A. López, R. Mart́ınez, J. A. Vera; Periodic orbits of the anisotropic Kepler problem
with quasihomogeneous potentials, International Journal of Bifurcation and Chaos 25 (2016),

no. 14, 1–7.



42 M. LI, Q. LIU EJDE-2021/63

[34] G. Maneff; La gravitation et le principe de l’égalité de l’action et de la réaction, Comptes

Rendus 178 (1924), 2159–2161.

[35] G. Maneff; Die gravitation und das prinzip von wirkung und gegenwirkung, Zeitschrift für
Physik volume 31 (1925), 786–802.

[36] G. Maneff; La gravitation et l’énergie au zéro, Comptes Rendus 190 (1930), 1374–1377.

[37] G. Maneff; Le principe de la moindre action et la gravitation, Comptes Rendus de I’Acadé
Bulgare des Sciences 190 (1930), 963–965.

[38] K. Meyer, G. Hall, O. Dan; Introduction to Hamiltonian dynamical systems and the n-body

problem, Springer New York, 2009.
[39] V. Mioc, C. Stoica; Discussion et résolution compléte du probléme des deux corps dans

le champ gravitationnel de maneff, Comptes Rendus de I’Acadé Bulgare des Sciences 320

(1995), 645–648.
[40] V. Mioc, C. Stoica; Discussion et résolution compléte du probléme des deux corps dans le

champ gravitationnel de maneff(ii), Comptes Rendus de I’Acadé Bulgare des Sciences 321
(1995), 961–964.

[41] M. Santoprete; Symmetric periodic solutions of the anisotropic Manev problem, Journal of

Mathematical Physics 43 (2002), no. 6, 3207–3219.
[42] F. Szenkovits, C. Stoica, V. Mioc; The Manev-Type problems: a topological view, Mathemat-

ica 2004 (2004), no. 1, 105–120.

[43] H. Tang, M. Xing, Q. Liu, G. Jiang; On the periodic orbits of four-particle time-dependent
fpu chains, Advances in Mathematical Physics 2018 (2018), 4526917.

[44] K. Tsetkova, V. Mioc; Manev’s field problemin contemporary science, American Institute of

Physics 1043 (2008), 137–141.
[45] F. Verhulst; Nonlinear differential equations and dynamical systems, World Publishing Corp.,

1993.

[46] C. Vidal; Periodic solutions for any planar symmetric perturbation of the Kepler problem,
Celestial Mechanics and Dynamical Astronomy 80 (2001), no. 2, 119–132.

Mengyuan Li

School of Mathematics and Computing Sciences, Guilin University of Electronic Tech-
nology, Guilin 541002, China.

Chengdu modern vocational & Technical school, Chengdu 610000, China

Email address: limy9595@163.com

Qihuai Liu

School of Mathematics and Computing Sciences, Guangxi Colleges and Universities
Key Laboratory of Data Analysis and Computation, Guilin University of Electronic

Technology, Guilin 541002, China.

Center for Applied Mathematics of Guangxi(Guangxi Normal University), Guilin 541004,
China

Email address: qhuailiu@gmail.com


	1. Introduction
	2. Applications related to spatial anisotropic problems
	2.1. The spatial anisotropic Manev problem
	2.2. The spatial anisotropic Kepler problem with perturbations

	3. Canonical transformations
	3.1. Spherical coordinates
	3.2. Spatial Delaunay elements

	4. Reduction of system on energy surface
	5. Standard form of averaging
	5.1. Calculation of integral for H10
	5.2. Calculation of integral for H11
	5.3. Averaged system

	6. Proof of main results
	Proof of Theorem ??
	6.1. Proof of Theorem ??

	7. Conclusion
	Acknowledgments

	References