Some Dynamic Aspects of a Sextic Galactic Potential in a Rotating Reference Frame

Abstract

This work aims to explore some dynamic aspects of the problem of star motion that is impacted by the rotation of the galaxy, which we model as a bisymmetric potential based on a two-dimensional harmonic oscillator with sextic perturbations. We demonstrate analytically that the motion is non-integrable when certain conditions are met. The analytical results for the non-integrability are confirmed by showing the irregularity of the behavior of the motion through utilizing the Poincaré surface of a section as a numerical method. The motion equilibrium positions are detected, and their stability is discussed. We show that the force generated by the rotating frame acts as a stabilizer for the maximum equilibrium points. We display graphically that the size of the stability regions relies on the angular velocity magnitude for the frame. Through the application of Lyapunov’s theorem, periodic solutions can be constructed which are close to the equilibrium positions. Furthermore, we demonstrate that there are one or two families of periodic solutions relying on whether the equilibrium point is a saddle or stable, respectively.

1. Introduction

Galactic dynamics, the study of the movements of galactic masses, is an important branch of astrophysics. It includes the investigations of stars, gases, and dark matter with the goal of understanding the main kinematics features of the three main types of galaxies (elliptical, spiral, and irregular) [1,2,3,4]. The theoretical bases for the subject were established approximately seven decades ago, when the study of dynamics was dominated by integrable and near-integrable systems. The study of galactic motion can be carried out globally. It can also be performed locally near the equilibrium points by using the perturbed oscillators. In the last few decades, several studies on galactic dynamics, mostly focused on elliptical galaxies, were published [5,6,7,8,9,10,11,12], and the existence of periodic solutions and their stability were investigated [13,14,15,16,17,18,19]. El-Sabaa et al. [20] also investigated the bifurcations of the Armbruster–Guckenheimer–Kim (AGK) potential by studying the topological type of the level sets of the integrable cases of the AGK galactic potential. This study was followed by Llibre et al. [21], who analyzed its global dynamics. Except for some studies on stellar orbits, the rotation of galaxies has been ignored in most studies of galactic dynamics [22]. This is in spite of the fact that an elliptical galaxy rotates with a slight angular velocity, and the force produced by the galaxy’s rotation will influence the dynamic (see, for example, [23,24]). This motivates the study of the rotation of a galaxy, modeled by a bisymmetric potential with different perturbations, and the effects of this rotation on the dynamics of star motion. Iñarrea et al. studied the stability of the motion equilibrium positions for a Hénon–Heiles system in a rotating reference frame [25]. Elmandouh in [26] studied the dynamic features of the AGK potential in a rotating reference frame. In that study, he examined the nonintegrability of the system to seek periodic solutions near the equilibrium points. Elmandouh et al. [27] investigated the integrability and stability of equilibria for a quartic galactic potential in a rotating reference frame. Lanchares et al. [28] applied Reeb’s theorem to prove the existence of periodic orbits in a rotating Hénon–Heiles system.

In the present work, we consider the plane motion rotation with a constant angular velocity $\mathrm{\Omega }$ of a galaxy approximately symmetrical around an axis normal to the plane of rotation. By reversing the normal axis if needed, we may assume that $\mathrm{\Omega }\ge 0$. This motion, similar to that of most elliptical galaxies, is described by a Hamiltonian system that involves a symmetric potential. The Hamiltonian we use is given by

$$H(x,y)=\frac{1}{2}(p_x^2+p_y^2)-\mathrm{\Omega }(x{p}_y-y{p}_x)+U(x,y),$$

(1)

where $p_x,p_y$ are the generalized momentum corresponding to the generalized coordinates $(x,y)$ and U is the four-parameter galactic potential function given by

$$U(x,y)=\frac{1}{2}(x^2+y^2)+\frac{1}{6}(a{x}^6+b{y}^6)+\frac{1}{2}{x}^2{y}^2(c{x}^2+d{y}^2),$$

(2)

while $a,b,c$, and d are free constants. Recently, Elsabaa et al. [29] studied the existence of periodic solutions and their stability by applying the averaging theorem for the sextic galactic potential while ignoring the rotation of the galaxy. This motivated us to study the influence of the galaxy rotations on some dynamics of the motion which were not considered before in the literature, such as integrability, the stability of the equilibrium positions, and the existence of periodic solutions near the equilibrium positions. The Hamilton canonical equations are

$$\begin{array}{cc}\hfill \dot{x}& =p_x+\mathrm{\Omega }y,\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \dot{y}& =p_y-\mathrm{\Omega }x,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\dot{p}}_x& =\mathrm{\Omega }{p}_y-x\left[1+a{x}^4+2c{x}^2{y}^2+d{y}^4\right],\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\dot{p}}_y& =-\mathrm{\Omega }{p}_x-y\left[1+b{y}^4+2d{x}^2{y}^2+c{x}^4\right].\hfill \end{array}$$

(3d)

As usual, the dots indicate that the derivative expressed is with respect to time. The equation of the motion (Equation (3)) has a Jacobi integral in the form

$$J_1:\frac{1}{2}(p_x^2+p_y^2)+\frac{1}{2}(x^2+y^2)+\frac{1}{6}(a{x}^6+b{y}^6)+\frac{1}{2}{x}^2{y}^2(c{x}^2+d{y}^2)=E,$$

(4)

where E is an arbitrary constant.

Additionally, the Hamiltonian in Equation (1) physically describes the motion of a particle under the combined action of the potential forces (velocity independent) obtained from the potential in Equation (2) and the gyroscopic forces (velocity-dependent) resulting from the rotation of the reference frame, giving the term $\mathrm{\Omega }(y{p}_x-x{p}_y)$.

This work is organized as follows. Section 2 contains the proof of the non-integrability of the Hamiltonian in Equation (1) and shows the behavior of the motion for two different values, leading to regular and irregular behavior of the motion. The equilibrium positions are determined in Section 3. In Section 4, we examine the stability of these equilibrium positions. We show the effect of the rotating reference frame on the stability and linear stability. The periodic solutions near the equilibrium positions are constructed, and the number of periodic solutions is linked to the stability of the equilibrium positions in Section 5.

2. Integrability Analysis

A 2D Hamiltonian system is integrable if it possesses two functionally independent integrals of motion $J_1$ and $J_2$ that are in involution (i.e., $\{J_1,J_2\}=0$), where $\{.,.\}$ represents Poisson brackets [30]. Integrable systems exhibit regularity, enabling the prediction of their behavior for all times. Their explicit solution can be found using quadratures, and they can be used to study non-integrable systems through perturbation theories [31]. The Hamiltonian system in Equation (1) is integrable if it has a complementary integral independent of the Jacobi integral in Equation (4). The problem of finding this integral is complicated due to the fact that its functional form is unknown. Details on constructing integrals of motion are mentioned in [32,33,34,35,36]. On the other hand, there are many methods that prove non-integrability, such as the Yoshida and Morales–Ramis criteria, which are based on the investigation of normal variation equations along a particular solution, and the Ziglin approach, which depends on the examination of the monodromy group associated with a variational equation [37,38,39]. Neither of these techniques work with Hamiltonian systems with linear momentum terms. As a consequence, we rewrote the Hamiltonian in Equation (1) appropriately.

We use the transformation

$$t\to {\alpha }^{-2}t,(x,y)\to (\alpha x,\alpha y),(p_x,p_y)\to ({\alpha }^3{p}_x,{\alpha }^3{p}_y).$$

(5)

Then, the Hamiltonian equations (Equation (3)) are transformed into

$$\begin{array}{cc}\hfill \dot{x}& =p_x+\frac{\mathrm{\Omega }}{{\alpha }^2}y,\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill \dot{y}& =p_y-\frac{\mathrm{\Omega }}{{\alpha }^2}x,\hfill \end{array}$$

(6b)

$$\begin{array}{cc}\hfill {\dot{p}}_x& =\frac{\mathrm{\Omega }}{{\alpha }^2}{p}_y-\frac{1}{{\alpha }^4}x-x[a{x}^4+2c{x}^2{y}^2+d{y}^4],\hfill \end{array}$$

(6c)

$$\begin{array}{cc}\hfill {\dot{p}}_y& =-\frac{\mathrm{\Omega }}{{\alpha }^2}{p}_x-\frac{1}{{\alpha }^4}y-y[b{y}^4+c{x}^4+2d{x}^2{y}^2].\hfill \end{array}$$

(6d)

When $\alpha$ tends toward infinity, the non-similarity system in Equation (6) becomes

$$\begin{array}{cc}\hfill \dot{x}& =p_x,\dot{y}=p_y,\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill {\dot{p}}_x& =-x[a{x}^4+2c{x}^2{y}^2+d{y}^4],\hfill \end{array}$$

(7b)

$$\begin{array}{cc}\hfill {\dot{p}}_y& =-y[b{y}^4+c{x}^4+2d{x}^2{y}^2].\hfill \end{array}$$

(7c)

The system in Equation (7) is a canonical Hamilton equation derived from the Hamiltonian

$$H_0=\frac{1}{2}(p_x^2+p_y^2)+\frac{1}{6}(a{x}^6+b{y}^6)+\frac{1}{2}{x}^2{y}^2(c{x}^2+d{y}^2).$$

(8)

We will study the non-integrability of the system in Equation (8) and show that the non-integrability of the Hamiltonian $H_0$ implies the non-integrability of the Hamiltonian H in Equation (1):

Theorem 1.

If the Hamiltonian $H_0$ is non-integrable, then the Hamiltonian H is non-integrable.

Proof. 

If the Hamiltonian H is integrable, then it has a complementary integral $J_2=J_2(x,y,p_x,p_y)$ independent of the Jacobi integral in Equation (4). By applying the transformation in Equation (5) to the complementary integral $J_2$ and taking the limit as $\alpha$ tending toward infinity, we obtain a complementary integral for $H_0$, confirming the integrability of $H_0$. The non-integrability of H follows by contradiction. □

Yoshida introduced a non-integrability theorem based on Ziglin’s analysis [37]. This theorem deals only with homogeneous potentials:

Theorem 2.

(Yoshida’s theorem) Let $H_0=\frac{1}{2}(p_x^2+p_y^2)+V_0(x,y)$ be a Hamiltonian system with homogeneous potential $V_0$ of a degree $k\ge 3$. Then, it is non-integrable if

$$\sigma \in ]-\infty ,0[\underset{j\in \mathbb{N}}{\cup }]j(j-1)\frac{k}{2}+j,j(j+1)\frac{k}{2}-j[,$$

(9)

where

$$\sigma =Trace\left({\mathit{Hess}}_{V_0}\left(\mathit{c}\right)\right)-k+1,$$

(10)

is the non-integrability coefficient and $\mathit{c}$ is the solution of the algebraic equation $\nabla V_0\left(\mathit{c}\right)=\mathit{c}$.

We are going to apply Theorem 2 to investigate the non-integrability of the Hamiltonian $H_0$, which is given by Equation (8). The potential $V_0(x,y)$ is

$$V_0(x,y)=\frac{1}{6}(a{x}^6+b{y}^6)+\frac{1}{2}{x}^2{y}^2(c{x}^2+d{y}^2).$$

(11)

The solutions to the algebraic equation $\nabla V_0(c_1,c_2)=(c_1,c_2)$ are as follows:

(1)

$c_1=0,c_2=\frac{1}{\sqrt[4]{b}}$, where $b>0$;

(2)

$c_1=\frac{1}{\sqrt[4]{a}},c_2=0$, where $a>0$.

For the first solution, the non-integrability coefficient is $\sigma =\frac{d}{b}$, while the non-integrability coefficient corresponding to the second solution is $\sigma =\frac{c}{a}$. Taking into account Yoshida’s theorem (Theorem 2), the Hamiltonian in Equation (8) is non-integrable if at least one of the following conditions is satisfied:

$$\begin{array}{cc}\hfill & \frac{d}{b}\in ]-\infty ,0[\cup ]j(3j-2),j(3j+2)[,j\in \mathbb{N},\hfill \\ \hfill & \frac{c}{a}\in ]-\infty ,0[\cup ]j(3j-2),j(3j+2)[,j\in \mathbb{N}.\hfill \end{array}$$

(12)

The next Theorem follows from Theorem 1 and Theorem 2:

Theorem 3. 

The galactic Hamiltonian in Equation (1) with the sextic galactic potential in Equation (2) is non-integrable if at least one of the conditions in Equation (12) is satisfied.

Remark 1. 

Theorem 3 provides the sufficient conditions for non-integrability. This means that if one of the conditions in Equation (2) is fulfilled, then the Hamiltonian in Equation (1) is non-integrable, and no further investigation is required. On the other hand, if neither conditions in Equation (2) hold, then the Hamiltonian in Equation (1) may be either integrable or non-integrable. In such a case, verifying integrability requires finding an additional integral functionally independent of the Jacobi integral in Equation (4).

It is worth motioning that Elmandouh introduced a theorem [34] based on studying the solvability of the differential Galois group associated with the normal variational equation, which can be applied effectively to give the necessary conditions for integrability and sufficient conditions for the non-integrability for a wide class of non-homogeneous potentials, including the galactic potential in Equation (2) as a special case.

Special cases:

(1)

If $a=b=c=d$, then the Hamiltonian in Equation (1) takes the form below in polar coordinates $(r,\theta )$:

$$H=\frac{1}{2}(p_r^2+\frac{1}{r^2}{p}_{\theta }^2)-\mathrm{\Omega }{p}_{\theta }+\frac{a}{6}{r}^6+\frac{r^2}{2}.$$

(13)

The Hamiltonian in Equation (13) is integrable since it has two independent first integrals of the motion. The first is the momentum $p_{\theta }=\alpha$ corresponding to the cyclic coordinate $\theta$, with the second being the Jacobi integral, which takes the form

$$\frac{1}{2}{p}_r^2+\frac{{\alpha }^2}{2r^2}+\frac{a}{6}{r}^6+\frac{r^2}{2}=E.$$

(14)

It should be noted that the constants $a,b,c$, and d in the present case do not satisfy Theorem 3 since $\frac{d}{b}=\frac{c}{a}=1$ does not fall in the non-integrability regions.

(2)

If $c=d=0$ and $\mathrm{\Omega }=0$, then the Hamiltonian in Equation (1) is separable into Cartesian coordinates. As a result of the separability of the Hamiltonian, the additional integral is quadratic in momentum. The integrals of the motion for this case are

$$\begin{array}{cc}\hfill & \frac{1}{2}{p}_x^2+\frac{1}{2}{x}^2+\frac{a}{6}{x}^6=\gamma ,\hfill \\ \hfill & \frac{1}{2}{p}_y^2+\frac{1}{2}{y}^2+\frac{b}{6}{y}^6=E-\gamma ,\hfill \end{array}$$

(15)

where $\gamma$ is a constant. It should be noticed that the non-integrability coefficient $\sigma$ is zero, and it does not fall within the non-integrability regions.

We present the Poincaré surface of a section for the Hamiltonian in Equation (1) with the potential function in Equation (2) when the parameters $a,b,c$, d, and $\mathrm{\Omega }$ take two different values. One of them satisfies Theorem 3, and the other does not agree with it. When $a=b=c=d$, and $\mathrm{\Omega }$ is arbitrary, the non-integrability coefficient $\sigma =1$ is outside the non-integrability region. The Poincaré surface of the section, which is illustrated by Figure 1a, seems to be regular and gives some indications about the integrability of the problem. This indication is confirmed by finding the complete set of the first integrals of the motion, in which one of them is a cyclic integral in polar coordinates $(r,\theta )$. If $a=1,b=1,c=5,d=7$, and $\mathrm{\Omega }$ is free, then Theorem (3) is met, and consequently, the Hamiltonian in Equation (1) is non-integrable. The Poincaré surface of the section corresponding to this case is irregular, which is usual for the non-integrability of the Hamiltonian in Equation (1). Hence, both the numerical and analytical results agree, and this shows the correctness of the obtained analytical results.

3. Equilibria

The equilibria of the system in Equation (3) can be obtained by setting the left-hand sides (the derivatives) in the Hamiltonian equations equal to zero to obtain

$$\begin{array}{cc}\hfill p_x& =-\mathrm{\Omega }y,p_y=\mathrm{\Omega }x,\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill \mathrm{\Omega }{p}_y& =x\left[1+a{x}^4+2c{x}^2{y}^2+d{y}^4\right],\hfill \end{array}$$

(16b)

$$\begin{array}{cc}\hfill \mathrm{\Omega }{p}_x& =-y\left[1+c{x}^4+2d{x}^2{y}^2+b{y}^4\right].\hfill \end{array}$$

(16c)

Substituting the two expressions in Equation (16a) into the Equations (16b) and (16c) gives

$$\begin{array}{cc}\hfill & x[a{x}^4+2c{x}^2{y}^2+d{y}^4+1-{\mathrm{\Omega }}^2]=0,\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill & y[b{y}^4+2d{x}^2{y}^2+c{x}^4+1-{\mathrm{\Omega }}^2]=0.\hfill \end{array}$$

(17b)

Rather than referring to the equilibrium positions as $(x,y,p_x,p_y)$, we indicate them as $(x,y)$, since $p_x,p_y$ can be directly evaluated by Equation (16a). By solving Equations (17a) and (17b), the equilibria are obtained in the next theorem:

Theorem 4. 

The Hamiltonian system in Equation (3), corresponding to the Hamiltonian in Equation (1), has the following as equilibrium positions:

1. 

If $a,b,c,d,1-{\mathrm{\Omega }}^2>0$, or if $a,b,c,d,1-{\mathrm{\Omega }}^2<0$, then ${\epsilon }_0=(0,0)$ is the only equilibrium point;

2. 

If $b({\mathrm{\Omega }}^2-1)>0$, then there is an equilibrium point ${\epsilon }_{1,2}=\left(0,\pm \sqrt[4]{\frac{{\mathrm{\Omega }}^2-1}{b}}\right)$;

3. 

If $a({\mathrm{\Omega }}^2-1)>0$, then there is an equilibrium point ${\epsilon }_{3,4}=\left(\pm \sqrt[4]{\frac{{\mathrm{\Omega }}^2-1}{a}},0\right)$;

4. 

If $b=a\ne c=d$ such that $(a+3c)({\mathrm{\Omega }}^2-1)>0$, then there is an equilibrium point ${\epsilon }_{5,6}=\left(\pm \sqrt[4]{\frac{{\mathrm{\Omega }}^2-1}{a+3c}},\pm \sqrt[4]{\frac{{\mathrm{\Omega }}^2-1}{a+3c}}\right)$.

It should be mentioned that there are some limiting cases in which either a, b, or $a+3c$ approaches zero, taking into account the equilibrium position existence condition. For example, if $b({\mathrm{\Omega }}^2-1)>0$, then ${\epsilon }_1$ blows up to infinity through the positive Y axis, while ${\epsilon }_2$ blows up to infinity through the negative Y axis. The same applies to the equilibrium positions ${\epsilon }_{3,4}$ via the positive and negative X axes. If the conditions in item 4 are satisfied, then equilibrium positions ${\epsilon }_{5,6}$ blow up to infinity through the line $y=x$ if the coordinates of the equilibrium point have the same sign or through the line $y=-x$ if the coordinates have a distinct sign. Moreover, if $a=b=c=d=0$, then the analysis of the dynamics reduces to that of a linear system which is readily understood. Dynamic analysis is more subtle if the constants are not zero simultaneously, as the nonlinear terms alter the behavior of the system.

4. Stability Analysis

First, we apply the Lagrange theorem [40], which has the advantage of yielding some information about the trapping and escape dynamics. To be a self-contained article, we recall the Lagrange theorem:

Theorem 5. 

(Lagrange theorem) [40] For a conservative system, if the potential energy has a strict minimum at some positions, then these positions are the positions of stable equilibrium points.

4.1. Application of the Lagrange Theorem

To apply the Lagrange Theorem 5, it is essential to find the effective potential of the problem. This requires a re-expression of the Hamiltonian in Equation (1) in terms of the generalized velocities $\dot{x},\dot{y}$ rather than the generalized momentum. Using the equations in Equation (3), we obtain

$$H(x,y)=\frac{1}{2}({\dot{x}}^2+{\dot{y}}^2)+V(x,y),$$

(18)

where $V(x,y)$ is the effective potential given by

$$V(x,y)=\frac{1}{2}(1-{\mathrm{\Omega }}^2)(x^2+y^2)+\frac{1}{6}(a{x}^6+b{y}^6)+\frac{1}{2}{x}^2{y}^2(c{x}^2+d{y}^2).$$

(19)

Notice that when $\mathrm{\Omega }=0$, the effective potential and potential in Equation (2) are identical. Since the system of equations $\frac{\partial V}{\partial x}=\frac{\partial V}{\partial y}=0$ is the same as in Equation (17), the equilibrium positions of the Hamiltonian system in Equation (1) are the critical points of the effective potential (Equation (19)). The nature of these positions is determined in Theorem 6 below by examining the Hessian matrix for the effective potential and using the Lagrange theorem. The following remark is needed:

Remark 2. 

If A is a non-singular $2\times 2$ matrix, then the matrix A is positive definite if and only if $det\left(A\right)>0, tr\left(A\right)>0$. The matrix A is negative definite if and only if $det\left(A\right)>0,tr\left(A\right)<0$.

Theorem 6. 

Under the existence conditions for the equilibrium points in Theorem 4, the following are satisfied:

1. 

The equilibrium ${\epsilon }_0$ is stable if $0\le \mathrm{\Omega }<1$. If $\mathrm{\Omega }>1$, then ${\epsilon }_0$ is an unstable maximum critical point;

2. 

The equilibrium ${\epsilon }_{1,2}$ is unstable when $\frac{d}{b}<1$. If $\frac{d}{b},\mathrm{\Omega }>1$, then ${\epsilon }_{1,2}$ is a stable position. If $\frac{d}{b}>1$ and $0\le \mathrm{\Omega }<1$, then the point ${\epsilon }_{1,2}$ is unstable and maximally critical;

3. 

The equilibrium ${\epsilon }_{3,4}$ is unstable when $\frac{c}{a}<1$. If $\frac{c}{a},\mathrm{\Omega }>1$, then ${\epsilon }_{3,4}$ is a stable position. The positions ${\epsilon }_{3,4}$ are maximally unstable if $\frac{c}{a}>1$ and $0\le \mathrm{\Omega }<1$;

4. 

The equilibrium ${\epsilon }_{5,6}$ is unstable when $(a-c)({\mathrm{\Omega }}^2-1)<0$. The case where $(a-c)({\mathrm{\Omega }}^2-1)>0$ implies that the point ${\epsilon }_{5,6}$ is stable if $a+c>0$. The positions ${\epsilon }_{5,6}$ are unstable and maximally critical if $(a-c)({\mathrm{\Omega }}^2-1)>0$ and $a+c<0$.

Proof. 

The Hessian matrix for the effective potential in Equation (19) takes the form

$$Hes{s}_V(x,y)=\left(\begin{array}{cc}1-{\mathrm{\Omega }}^2+5a{x}^4+6c{x}^2{y}^2+d{y}^4& 4xy(c{x}^2+d{y}^2)\\ \\ 4xy(c{x}^2+d{y}^2)& 1-{\mathrm{\Omega }}^2+5b{y}^4+6d{x}^2{y}^2+c{x}^4\end{array}\right).$$

(20)

We examine the Hessian matrix at the equilibrium positions calculated in Theorem 4. For ${\epsilon }_0=(0,0)$, the Hessian matrix in Equation (20) takes the form

$$\begin{array}{c}\hfill Hes{s}_V\left({\epsilon }_0\right)=\left(\begin{array}{cc}1-{\mathrm{\Omega }}^2& 0\\ 0& 1-{\mathrm{\Omega }}^2\end{array}\right),\end{array}$$

(21)

with $det\left(Hes{s}_V\left({\epsilon }_0\right)\right)={(1-{\mathrm{\Omega }}^2)}^2>0$. Using Remark 2, the point ${\epsilon }_0$ is either a local minimum if $0\le \mathrm{\Omega }<1$ or local maximum if $\mathrm{\Omega }>1$. Based on the Lagrange theorem, ${\epsilon }_0$ is stable if $0\le \mathrm{\Omega }<1$ and unstable if $\mathrm{\Omega }>1$. The Hessian matrix in Equation (20) at ${\epsilon }_{1,2}$ is

$$Hes{s}_V\left({\epsilon }_{1,2}\right)=\left(\begin{array}{cc}({\mathrm{\Omega }}^2-1)(\frac{d}{b}-1)& 0\\ 0& 4({\mathrm{\Omega }}^2-1)\end{array}\right),$$

(22)

whose determinant is $4{({\mathrm{\Omega }}^2-1)}^2(\frac{d}{b}-1)$ and whose trace is $({\mathrm{\Omega }}^2-1)(\frac{d}{b}+3)$. This implies the following:

• If $\frac{d}{b}<1$, then $Hes{s}_V\left({\epsilon }_{1,2}\right)$ is an indefinite matrix, and the point ${\epsilon }_{1,2}$ is a saddle and hence unstable.
• If $\frac{d}{b}>1$, the point ${\epsilon }_{1,2}$ is a local minimum and hence stable if $\mathrm{\Omega }>1$. If $0\le \mathrm{\Omega }<1$, then it is a local maximum and hence unstable.

Similarly, the case for ${\epsilon }_{3,4}$ with $a,c$ replacing $b,d$. For the equilibrium point ${\epsilon }_{5,6}$, the Hessian matrix (Equation (20)) is

$$Hes{s}_V\left({\epsilon }_{5,6}\right)=\frac{4({\mathrm{\Omega }}^2-1)}{a+3c}\left(\begin{array}{cc}(a+c)& \pm 2c\\ \pm 2c& (a+c)\end{array}\right),$$

(23)

where the off-diagonal entries have the same sign. The determinant of $Hes{s}_V\left({\epsilon }_{5,6}\right)$ is $4({\mathrm{\Omega }}^2-1)(a-c)$. If $({\mathrm{\Omega }}^2-1)(a-c)<0$, then the equilibrium point ${\epsilon }_{5,6}$ is a saddle and hence unstable. If $({\mathrm{\Omega }}^2-1)(a-c)>0$, then we check the sign of the trace which follows the sign of $(a+c)$. If $a+c>0$, then the point ${\epsilon }_{5,6}$ is a local minimum and hence stable. If $a+c<0$, then the case where $(a-c)$ and $({\mathrm{\Omega }}^2-1)$ are both negative gives ${\epsilon }_{5,6}$ as a local maximum and is thus unstable. The case where both $(a-c)$ and $({\mathrm{\Omega }}^2-1)$ are positive and $a+c<0$ leads to a contradiction, so it never happens. □

It should be noticed that in all the above cases, the Hessian matrix becomes null when $\mathrm{\Omega }=1$ and no information about the equilibrium positions can be derived. Consequently, we exclude the case where $\mathrm{\Omega }=1$ from further discussion, as it demands a special investigation outside the scope of this work.

Notice that the degree of instability for the maximum equilibrium points is two, and according to Lord Kelvin’s proposition [41], with the addition of gyroscopic forces, these positions remain stable when they are the minimum for the effective potential in Equation (19), while the unstable maximum equilibrium points for the effective potential can be stabilized by the forces resulting from the rotating frame (gyroscopic forces). Consequently, we must apply the linear approximation for the stability of these positions [40]:

4.2. Linear Stability

This section completes the study in Section 4.1 by investigating the linear stability of the equilibrium positions given in Theorem 4. The linearization system corresponding to the Hamiltonian system in Equation (3) takes the form

$$\dot{\mathit{z}}=\mathcal{J}\mathit{z},$$

(24)

where $\mathit{z}=(z_1,z_2,z_3,z_4)$ are the perturbed variables for $x,y,p_x$, and $p_y$, respectively. The matrix $\mathcal{J}$ is the Jacobian matrix evaluated at the equilibrium position $(x_0,y_0)$ given by

$$\mathcal{J}(x_0,y_0)=\left(\begin{array}{cccc}0& \mathrm{\Omega }& 1& 0\\ -\mathrm{\Omega }& 0& 0& 1\\ A& B& 0& \mathrm{\Omega }\\ B& C& -\mathrm{\Omega }& 0\end{array}\right).$$

(25)

where $A=-5a{x}_0^4-6c{x}_0^2{y}_0^2-d{y}_0^4-1,B=-4x_0{y}_0(c{x}_0^2+d{y}_0^2)$, and $C=-5b{y}_0^4-c{x}_0^4-6d{x}_0^2{y}_0^2-1$. The eigenvalues of the Jacobian matrix in Equation (25) are the roots of the characteristic equation

$${\lambda }^4+M{\lambda }^2+N=0,$$

(26)

where M and N are given by

$$M=2{\mathrm{\Omega }}^2-(A+C),N=({\mathrm{\Omega }}^2+A)({\mathrm{\Omega }}^2+C)-B^2.$$

(27)

The equilibrium position $(x_0,y_0)$ is linearly stable if the zeros of the characteristic equation (eigenvalues of the Jacobian matrix in Equation (25)) have pure imaginary roots, and the Jacobian matrix is semi-simple. Thus, the conditions for linear stability are

$$N\ge 0,M\ge 0,M^2-4N\ge 0.$$

(28)

It is worth mentioning that if at least one of the conditions in Equation (28) is not satisfied at $(x_0,y_0)$, then the equilibrium position $(x_0,y_0)$ is unstable, and the problem does not require further investigation. Based on the Lyapunov theorem (see, for example, [41]), the instability of the equilibrium position $(x_0,y_0)$ in a linear approximation is still valid when the nonlinear terms in the Hamilton system in Equation (3) are considered, while the stability of the equilibrium position $(x_0,y_0)$ is not necessarily still valid when the nonlinear terms are taken into account. This means the linear approximation analysis gives the necessary conditions for stability and sufficient conditions for instability.

According to Theorem 6, the equilibrium point ${\epsilon }_0$ is a maximum critical point for the effective potential in Equation (19) if $\mathrm{\Omega }>1$ and the eigenvalues of the Jacobian matrix are $\pm (\mathrm{\Omega }-1)i$ and $\pm (\mathrm{\Omega }+1)i$. Hence, ${\epsilon }_0$ is a linearly stable maximum. This proves the first statement in the following theorem:

Theorem 7. 

For the equilibrium positions of the Hamiltonian system (3), we have

1. 

If all the coefficients $a,b,c,d<0$ and $\mathrm{\Omega }>1$, then ${\epsilon }_0$ is a linearly stable maximum.

2. 

Assuming $b({\mathrm{\Omega }}^2-1)>0$, the positions ${\epsilon }_{1,2}$ are linearly stable maximums if $d\in [b{\mu }_1,b[$ and $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$. Otherwise, they are Lyapunov unstable.

3. 

Assuming $a({\mathrm{\Omega }}^2-1)>0$, the positions ${\epsilon }_{3,4}$ are linearly stable maximums if $c\in [a{\mu }_1,a[$ and $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$. Otherwise, they are Lyapunov unstable.

where ${\mu }_1=(5-{\mathrm{\Omega }}^2-8\mathrm{\Omega }\sqrt{1-{\mathrm{\Omega }}^2})/(1-{\mathrm{\Omega }}^2)$.

Proof. 

By Theorem 6, the equilibrium positions ${\epsilon }_{1,2}$ are the maximums for the effective potential in Equation (19) if $\frac{d}{b}>1,0\le \mathrm{\Omega }<1$. To maintain the equilibrium condition, since $\mathrm{\Omega }<1$, b must be negative, then $d<b$. Let $\mu =\frac{d}{b}$. By calculating the expressions for M and N at the equilibrium positions ${\epsilon }_{1,2}$, we have

$$\begin{array}{cc}\hfill M& =7{\mathrm{\Omega }}^2-3+\mu ({\mathrm{\Omega }}^2-1),\hfill \end{array}$$

(29a)

$$\begin{array}{cc}\hfill N& =4{({\mathrm{\Omega }}^2-1)}^2(\mu -1),\hfill \end{array}$$

(29b)

$$\begin{array}{cc}\hfill M^2-4N& ={\left[7{\mathrm{\Omega }}^2-3+\mu ({\mathrm{\Omega }}^2-1)\right]}^2-16{({\mathrm{\Omega }}^2-1)}^2(\mu -1).\hfill \end{array}$$

(29c)

The equilibrium positions ${\epsilon }_{1,2}$ are linearly stable if the conditions in Equation (28) are satisfied, taking into account the existence of the condition $b<0,0\le \mathrm{\Omega }<1$. The two conditions $M\ge 0$ and $N\ge 0$ are satisfied together if

$$1<\mu \le \frac{3-7{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}^2-1}.$$

(30)

To guarantee the correctness of the last inequality, we find the possible values of $\mathrm{\Omega }$ such that $1<\frac{3-7{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}^2-1}$ is verified. Thus, we have $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$. Let $f\left(\mu \right)=M^2-4N$ and $I=]1,\frac{3-7{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}^2-1}]$.

Direct calculations show that $\frac{df}{d\mu }<0$ for all $\mu \in I$, and

$$f\left(1\right)=16{(2{\mathrm{\Omega }}^2-1)}^2>0,f\left(\frac{3-7{\mathrm{\Omega }}^2}{{\mathrm{\Omega }}^2-1}\right)=64(2{\mathrm{\Omega }}^2-1)({\mathrm{\Omega }}^2-1)<0.$$

Thus, the polynomial $f\left(\mu \right)$ has a unique root ${\mu }_1=\frac{5-{\mathrm{\Omega }}^2-8\mathrm{\Omega }\sqrt{1-{\mathrm{\Omega }}^2}}{1-{\mathrm{\Omega }}^2}$ in the interval I. Hence, the equilibrium position ${\epsilon }_{1,2}$ is a linearly stable maximum if

$$d<b<0,\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[,1<\frac{d}{b}\le \frac{5-{\mathrm{\Omega }}^2-8\mathrm{\Omega }\sqrt{1-{\mathrm{\Omega }}^2}}{1-{\mathrm{\Omega }}^2}.$$

(31)

Analogous procedures can be applied to investigate the linear stability of the equilibrium position ${\epsilon }_{3,4}$. □

The regions of linear stability of the equilibrium point ${\epsilon }_{1,2}$ are delimited by the two straight lines $b=d$ and $d={\mu }_1\left(\mathrm{\Omega }\right)b$, where ${\mu }_1\left(\mathrm{\Omega }\right)$ is the slope of the second line as well as the region in which the equilibrium positions ${\epsilon }_{1,2}$ are the maximums, as determined by the two lines $d=0,b=d$. Hence, the size of the region of linear stability relies on the value of the angular velocity $\mathrm{\Omega }$. Figure 2 outlines the zone of linear stability in yellow, which grows as the value of the angular velocity $\mathrm{\Omega }$ increases, and it reaches the maximum size when $\mathrm{\Omega }\to 1$ (i.e., it covers all the regions of the maximum positions in gray for ${\epsilon }_{1,2}$). This means the force resulting from the rotating frame controls the regions of stability and instability. On another side, if at least one condition in Equation (31) is not fulfilled, then the equilibrium position ${\epsilon }_{1,2}$ becomes Lyapunov unstable. The same discussion and similar graphs can be constructed for the equilibrium point ${\epsilon }_{3,4}$ using $c,a$ instead of $d,b$.

According to Theorem 6, the equilibrium positions ${\epsilon }_{5,6}$ are the maximums for the effective potential in Equation (19) if $(a-c)({\mathrm{\Omega }}^2-1)>0$ and $a+c<0$. The following theorem shows the conditions for their linear stability, taking into account the existence condition $(a+3c)({\mathrm{\Omega }}^2-1)>0$.

Theorem 8. 

Assuming $(a+3c)({\mathrm{\Omega }}^2-1)>0$, the equilibrium positions ${\epsilon }_{5,6}$ of the Hamiltonian system in Equation (3) are linearly stable maximums if $\mathrm{\Omega }\in ]\sqrt{\frac{1}{2}},1[$ and $c\in ]a,{\mu }_{2}a]$, where

$${\mu }_2=\frac{8{\mathrm{\Omega }}^2-11{\mathrm{\Omega }}^4+8\mathrm{\Omega }(1-{\mathrm{\Omega }}^2)\sqrt{1-{\mathrm{\Omega }}^2}}{37{\mathrm{\Omega }}^4-44{\mathrm{\Omega }}^2+16}.$$

Proof. 

Considering the equilibrium condition $(a+3c)({\mathrm{\Omega }}^2-1)>0$, the quantities $a+3c$ and $({\mathrm{\Omega }}^2-1)$ are both positive or both negative. Therefore, the equilibrium positions ${\epsilon }_{5,6}$ are the maximums for the effective potential in Equation (19) if all the quantities $(a-c),(a+3c)$, and $a+c$ are negatives (the case where $(a-c)$ and $(a+3c)$ are both positive and $(a+c)$ is negative will lead to a contradiction and never happen). To verify the linear stability of ${\epsilon }_{5,6}$, we compute the corresponding M and N in Equation (27) to obtain

$$\begin{array}{cc}\hfill N& =\frac{16(a-c){({\mathrm{\Omega }}^2-1)}^2}{a+3c},\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill M& =\frac{4(a(3{\mathrm{\Omega }}^2-2)+c(5{\mathrm{\Omega }}^2-2))}{a+3c}=\frac{4(3a+5c)}{a+3c}{\mathrm{\Omega }}^2-\frac{8(a+c)}{a+3c},\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill M^2-4N& =\frac{16((37{\mathrm{\Omega }}^4-44{\mathrm{\Omega }}^2+16)c^2+2ac{\mathrm{\Omega }}^2(11{\mathrm{\Omega }}^2-8)+a^2{\mathrm{\Omega }}^2(5{\mathrm{\Omega }}^2-4))}{{(a+3c)}^2}.\hfill \end{array}$$

(32c)

Note that N is positive. To check the signs of the other quantities, let $\mu =\frac{c}{a}$. Since $a+3c<0,a-c<0$, then $a<0$, and it easily follows that $\mu \in ]-\frac{1}{3},1[$:

• To find the conditions for $M^2-4N\ge 0$, we write $M^2-4N$ as a function in the two variables $\mu$ and $\mathrm{\Omega }$ to obtain

  $$M^2-4N=\frac{16((37{\mathrm{\Omega }}^4-44{\mathrm{\Omega }}^2+16){\mu }^2+2\mu {\mathrm{\Omega }}^2(11{\mathrm{\Omega }}^2-8)+{\mathrm{\Omega }}^2(5{\mathrm{\Omega }}^2-4))}{{(1+3\mu )}^2}.$$
  For each fixed $\mathrm{\Omega }$, the sign of $M^2-4N$ is determined by the sign of

  $$f\left(\mu \right)=(37{\mathrm{\Omega }}^4-44{\mathrm{\Omega }}^2+16){\mu }^2+2{\mathrm{\Omega }}^2(11{\mathrm{\Omega }}^2-8)\mu +{\mathrm{\Omega }}^2(5{\mathrm{\Omega }}^2-4).$$
  Since the discriminate of $f\left(\mu \right)$ is $-256{\mathrm{\Omega }}^2{({\mathrm{\Omega }}^2-1)}^3>0$, there are two different roots for $f\left(\mu \right)$, given by

  $${\mu }_{1,2}=\frac{8{\mathrm{\Omega }}^2-11{\mathrm{\Omega }}^4\pm 8\mathrm{\Omega }(1-{\mathrm{\Omega }}^2)\sqrt{1-{\mathrm{\Omega }}^2}}{37{\mathrm{\Omega }}^4-44{\mathrm{\Omega }}^2+16}.$$
  It is easy to verify that $\frac{-1}{3}<{\mu }_1<{\mu }_2\le 1$ for $\mathrm{\Omega }\in [0,1[$, where ${\mu }_2=1$ at $\mathrm{\Omega }=\frac{1}{\sqrt{2}}$. For each fixed $\mathrm{\Omega }$ such that $\mathrm{\Omega }\ne \frac{1}{\sqrt{2}}$, the function is $f\left(\mu \right)$, and hence $M^2-4N$ is non-negative on $]\frac{-1}{3},{\mu }_1]\cup [{\mu }_2,1[$ and negative on $]{\mu }_1,{\mu }_2[$. Therefore, we have to exclude the interval $]{\mu }_1,{\mu }_2[$ from the possible values of $\mu$.
• On the interval $]\frac{-1}{3},{\mu }_1]\cup [{\mu }_2,1[$, we investigate the sign of M. By writing M as a function in the two variables $\mu$ and $\mathrm{\Omega }$, we obtain

  $$M=\frac{4(3+5\mu )}{1+3\mu }{\mathrm{\Omega }}^2-\frac{8(1+\mu )}{1+3\mu }.$$

  (33)
  For each $\mu \in ]-\frac{1}{3},1[$, all the quantities $1+\mu ,1+3\mu$, and $3+5\mu$ are positive, and the map $g\left(\mu \right)=\frac{2(1+\mu )}{3+5\mu }$ is decreasing on $\mu \in ]-\frac{1}{3},1[$. For a fixed $\mu$, the quantity $M\ge 0$ if and only if ${\mathrm{\Omega }}^2\ge \frac{2(1+\mu )}{3+5\mu }=g\left(\mu \right)$. Direct computations show the following:

  –

  At the fixed value ${\mu }_1$, we have ${\mathrm{\Omega }}^2<g\left({\mu }_1\right)$ for all $\mathrm{\Omega }\in [0,1[$. Therefore, for each $\mu \in ]-\frac{1}{3},{\mu }_1]$ and $\mathrm{\Omega }\in [0,1[$, we have ${\mathrm{\Omega }}^2<g\left({\mu }_1\right)\le g\left(\mu \right)$, and then $M<0$. Therefore, we must also exclude the interval $]\frac{-1}{3},{\mu }_1]$ from the possible values for $\mu$.

  –

  At the fixed value ${\mu }_2$, we have ${\mathrm{\Omega }}^2\ge g\left({\mu }_2\right)$ for each $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$. Therefore, for each $\mu \in [{\mu }_2,1[$ and $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$, we have ${\mathrm{\Omega }}^2\ge g\left({\mu }_2\right)\ge g\left(\mu \right)$, and then $M\ge 0$. The case where ${\mathrm{\Omega }}^2=\frac{1}{\sqrt{2}}$ yields ${\mu }_2=1$ and must be excluded.

  –

  At the fixed point $\mu =1$, we have ${\mathrm{\Omega }}^2<g\left(1\right)$ for $\mathrm{\Omega }\in ]0,\frac{1}{\sqrt{2}}[$. Therefore, for each $\mu \in ]{\mu }_2,1[$ and $\mathrm{\Omega }\in ]0,\frac{1}{\sqrt{2}}[$, we have ${\mathrm{\Omega }}^2<g\left(1\right)<g\left(\mu \right)$, and then $M<0$. Therefore, we must also exclude the interval $\mathrm{\Omega }\in ]0,\frac{1}{\sqrt{2}}[$ from the possible values of $\mathrm{\Omega }$. Note that even though $\mu =1$ is not included in the interval under consideration, $g\left(\mu \right)$ is continuous at $\mu =1$, so we can study the relation between $\mathrm{\Omega }$ and $\mu$ when $\mu$ tends toward 1.

  Therefore, all the three quantities $N,M,M^2-4M$ are non-negative if $\mu \in [{\mu }_2,1[,$ and $\mathrm{\Omega }\in ]\frac{1}{\sqrt{2}},1[$.

□

The zone of linear stability for the motion equilibrium positions is specified by the two lines $c=a$ and $c={\mu }_2\left(\mathrm{\Omega }\right)a$, displayed in yellow in Figure 3, as well as the zone in which these positions are the maximums, which are limited by the lines $c=-\frac{1}{3a}$ and $c=-a$, as clarified in yellow in Figure 3. The size of the linear stability relies on the value of the angular velocity $\mathrm{\Omega }$. As shown in Figure 3b,c, the region of linear stability grows as the angular velocity tends toward one, and at the same time, the size of the Lyapunov instability decreases. Moreover, when $\mathrm{\Omega }$ tends toward one, the linear stability covers the regions in which ${\epsilon }_{5,6}$ are the maximums. Hence, the force resulting from the rotating reference frame can be used as a stabilizer for the unstable maximum equilibrium points.

5. Periodic Motion

In this section, periodic solutions for the Hamiltonian system in Equation (3) around its equilibrium positions are discussed based on the Lyapunov theorem [42], which has been successfully applied in several works. Yehia [43] presented the periodic solutions near the equilibrium position for mechanical systems in general, while El-Sabaa introduced periodic solutions for the motion of a rigid body’s dynamics in a central Newtonian field [44]. Elmandouh [26] found periodic solutions for the AGK potential nearest to its equilibrium points. El-Sabaa et al. [29] constructed periodic solutions for the sextic galactic potential in a non-rotating reference frame.

The Lagrangian equations of the motions are

$$\ddot{x}-2\mathrm{\Omega }\dot{y}+V_x=0,\ddot{y}+2\mathrm{\Omega }\dot{x}+V_y=0,$$

(34)

where V is given by Equation (19). The equations of the motion in Equation (34) has a first integral

$$\frac{1}{2}({\dot{x}}^2+{\dot{y}}^2)+V=E,$$

(35)

where E is an arbitrary constant. This integral is called the Jacobi integral. Since the equilibrium point $(x_0,y_0)$ is a critical point for the effective potential V, we have

$$V_x(x_0,y_0)=V_y(x_0,y_0)=0,E_0=V(x_0,y_0),$$

(36)

where $E_0$ denotes the Jacobi integral value at the point $(x_0,y_0)$. To describe the perturbed motion about $(x_0,y_0)$, we introduce

$$x=x_0+p,y=y_0+q,$$

(37)

where p and q are perturbation variables. Notice that the expressions in Equation (37) transform the equilibrium point $(x_0,y_0)$ to the origin of the coordinates. When inserting the expressions in Equation (37) into the equations of motion (34), we have

$$\ddot{p}-2\mathrm{\Omega }\dot{q}+V_{xx}(x_0,y_0)p+V_{xy}(x_0,y_0)q=0,\ddot{q}+2\mathrm{\Omega }\dot{p}+V_{xy}(x_0,y_0)p+V_{yy}(x_0,y_0)q=0.$$

(38)

Using a new variable $\zeta$ instead of the time through the time transformation yields

$$\zeta =\gamma t,$$

(39)

where $\gamma$ is a free parameter. By applying the transformation in Equation (39), the system in Equation (38) is converted into

$$\begin{array}{cc}\hfill & {\gamma }^2{p}^{″}-2\gamma \mathrm{\Omega }{q}^{\prime }+V_{xx}(x_0,y_0)p+V_{xy}(x_0,y_0)q=0,\hfill \\ \hfill & {\gamma }^2{q}^{″}+2\gamma \mathrm{\Omega }{p}^{\prime }+V_{xy}(x_0,y_0)p+V_{yy}(x_0,y_0)q=0,\hfill \end{array}$$

(40)

where ${}^{\prime }$ refers to the derivative with respect to $\zeta$. Following the Lyapunov theorem [42], we assume the periodic solutions as a power series of n which, in our case, depends on the Jacobi integral value (i.e., it depends on E). Assume that

$$p=\sum _{j=1}^{\infty }{n}^j{x}_j,q=\sum _{j=1}^{\infty }{n}^j{y}_j,E=E_0+\sum _{j=2}^{\infty }{n}^j{E}_j,$$

(41)

where $E_i$ are constants and $x_i$ and $y_i$ are $T-$ periodic functions of t with a period $\gamma$ that can be written as

$$T=\frac{2\pi }{\gamma }=\frac{2\pi }{{\gamma }_0}\left(1+\sum _{j=2}^{\infty }{n}^j{T}_j\right).$$

(42)

By substituting the expressions in Equations (41) and (42) into Equation (40) and finding the first order approximation, we have

$$\begin{array}{cc}\hfill & {\gamma }_0^2{x}_1^{″}+2{\gamma }_0\mathrm{\Omega }{y}_1^{\prime }+V_{xy}(x_0,y_0)y_1+V_{yy}(x_0,y_0)x_1=0,\hfill \\ \hfill & {\gamma }_0^2{y}_1^{″}-2\mathrm{\Omega }{\gamma }_0{x}_1^{\prime }+V_{xx}(x_0,y_0)y_1+V_{xy}(x_0,y_0)x_1=0.\hfill \end{array}$$

(43)

Our goal is accomplished by looking for solutions to Equation (43) in the form of

$$\begin{array}{c}{x}_1=a_0+a_{1}cos\zeta +a_{2}sin\zeta ,\hfill \\ y_1=b_0+b_{1}cos\zeta +b_{2}sin\zeta ,\hfill \end{array}$$

(44)

where $a_i,b_i,i=0,1,2$ are the free constants that need to be determined. When putting the expressions in Equation (44) into Equation (43), we obtain

$$\begin{array}{cc}\hfill x\left(\zeta \right)& =x_0+n\left[k_1{V}_{xy}(x_0,y_0)cos\zeta +(({\gamma }_0^2-V_{xx}(x_0,y_0))k_1-2\mathrm{\Omega }{\gamma }_0{k}_2)sin\zeta \right],\hfill \\ \hfill y\left(\zeta \right)& =y_0+n\left[k_1(({\gamma }_0^2-V_{yy}(x_0,y_0))-2k_2\mathrm{\Omega }{\gamma }_0)cos\zeta +k_2{V}_{xy}(x_0,y_0)sin\zeta \right],\hfill \\ \hfill E& =E_0+n^2\left[\frac{k_1^2+k_2^2}{4}{\gamma }_0^6-2\mathrm{\Omega }{k}_1{k}_2{\gamma }_0^5+\frac{{\gamma }_0^4}{4}\left(k_1^2\right(4{\mathrm{\Omega }}^2+V_{xx}(x_0,y_0)-2V_{yy}(x_0,y_0)\right.\hfill \\ \hfill & +k_2^2(4{\mathrm{\Omega }}^2-2V_{xx}(x_0,y_0)+2V_{yy}(x_0,y_0))-\frac{1}{4}{\gamma }_0^2(-k_2^2{V}_{xx}(x_0,y_0)+2((k_1^2+k_2^2)V_{yy}(x_0,y_0)\hfill \\ \hfill & -2k_2^2{\mathrm{\Omega }}^2)V_{xx}-k_1^2{V}_{yy}^2(x_0,y_0)-4k_1^2{\mathrm{\Omega }}^2{V}_{yy}(x_0,y_0)-3(k_1^2+k_2^2)V_{xy}^2(x_0,y_0)\left)\right)\hfill \\ \hfill & \left.+\frac{(V_{xx}(x_0,y_0)V_{yy}^2(x_0,y_0)}{4}(k_1^2{V}_{yy}(x_0,y_0)+k_2^2{V}_{xx}(x_0,y_0)\right],\hfill \end{array}$$

(45)

where $k_1$ and $k_2$ are arbitrary constants introduced instead of the old ones for the suitability, while ${\gamma }_0$ is the frequency, and it takes the value

$${\gamma }_0=\sqrt{\frac{V_{yy}(x_0,y_0)+V_{xx}(x_0,y_0)+4{\mathrm{\Omega }}^2\pm \sqrt{{(V_{yy}(x_0,y_0)+V_{xx}(x_0,y_0)+4{\mathrm{\Omega }}^2)}^2-4\mathcal{R}}}{2}},$$

(46)

where $\mathcal{R}$ is given by

$$\mathcal{R}=V_{xx}(x_0,y_0)V_{yy}(x_0,y_0)-V_{xy}^2(x_0,y_0).$$

(47)

Based on the values $V_{xx}(x_0,y_0),V_{yy}(x_0,y_0),V_{xy}(x_0,y_0)$, and $\mathrm{\Omega }$, there are single or double values of the frequency ${\gamma }_0$. Let us clarify this point:

• If $\mathcal{R}=V_{xx}(x_0,y_0)V_{yy}(x_0,y_0)-V_{xy}^2(x_0,y_0)>0$, this means $(x_0,y_0)$ is an extreme point for the effective potential in Equation (19). Thus, we have two frequencies if the condition

  $${\mathrm{\Omega }}^2>\frac{1}{4}[2\sqrt{\mathcal{R}}-V_{yy}(x_0,y_0)-V_{xx}(x_0,y_0)],$$

  (48)

  is verified. This condition is always satisfied when $(x_0,y_0)$ is the minimum value of the effective potential. Thus, there are two values for the frequency. Consequently, there are two periodic families of orbits depending on the parameter n about the stable minimum point of the effective potential.
• If $\mathcal{R}=V_{xx}(x_0,y_0)V_{yy}(x_0,y_0)-V_{xy}^2(x_0,y_0)<0$, then the point $(x_0,y_0)$ is a saddle point for the effective potential in Equation (19). There is a unique value for the frequency ${\gamma }_0$ which is obtained from Equation (46), which takes the positive sign. Hence, there is exactly one periodic orbit around the saddle point for the equilibrium point.

We outlined the obtained results graphically by introducing certain values for the included parameters. The equilibrium point ${\epsilon }_0=(0,0)$ is a stable minimum critical point for the effective potential in Equation (19) if we select $\mathrm{\Omega }=0.5,n={10}^{-5},k_1=0.5$, and $k_2=1.5$. For these values, Equation (46) gives two real frequencies: ${\gamma }_0=1.5$ and ${\gamma }_0=0.5$. Thus, there are two periodic orbits close to the equilibrium point ${\epsilon }_0$, as outlined in Figure 4a. For the equilibrium point ${\epsilon }_{1,2}$, we choose $\mathrm{\Omega }=1.5,b=0.2$, which guarantees the existence of this equilibrium alongside $n={10}^{-5},k_1=0.5$. There are two possible choices for the parameter d:

• If $d=0.1$, then the equilibrium point ${\epsilon }_{1,2}$ is a saddle point for the effective potential since $\frac{d}{b}=0.5<1$. Equation (46) gives a single frequency ${\gamma }_0=3.688455$. Thus, one periodic orbit is close to ${\epsilon }_{1,2}$, as outlined in Figure 4b.
• If $d=0.4$, then $\frac{d}{b}=2>1$, and so the equilibrium point ${\epsilon }_{1,2}$ is a stable minimum critical point for the effective potential in Equation (19). Equation (46) provides two real values of the frequency: ${\gamma }_0=3.85078$ and ${\gamma }_0=0.649219$. Hence, there are two periodic families corresponding to each frequency, as illustrated by Figure 4c.

Figure 4. Periodic solutions close to the equilibrium points ${\epsilon }_0$ and ${\epsilon }_{1,2}$. The green solid circle indicates the equilibrium position.

Figure 4. Periodic solutions close to the equilibrium points ${\epsilon }_0$ and ${\epsilon }_{1,2}$. The green solid circle indicates the equilibrium position.

6. Discussion

This work explored some aspects of motion under the sextic galactic potential. More specifically, we studied the influence of the forces resulting from a rotating reference frame on the dynamics. The non-integrability of the problem was examined using Yoshida’s analysis. This analysis applies only to the homogeneous potential, so we applied certain similar transformations to convert the original non-homogeneous potential to a homogeneous one. We also used a Poincare surface section as a numerical tool to confirm the obtained analytical results for two values of the parameters: one satisfying the non-integrability Theorem 3 and the other not satisfying it. The stability of the equilibrium position was investigated by applying the Lagrange theorem (Theorem 5), which implies the stability of the minimum positions and the instability of the maximum positions of the effective potentials. As a result, the degree of instability of the maximum equilibrium positions was two. The gyroscopic forces, which resemble the forces resulting from the rotating reference frame, can be utilized as a stabilizer for these unstable positions. Therefore, we applied the linear approximation method, which was also employed to illustrate the dependence of the linear stability regions on the angular velocity. Some figures in the plane of the parameters were used to describe the dependence of the stability zones on the angular velocity. The Lyapunov theorem has been applied to construct periodic solutions that are almost close to the equilibrium positions. We also proved that near the stable equilibrium positions, there are two periodic families as well as one periodic family of orbits near the unstable equilibrium positions.