Perturbed Planar Restricted Four-Body Problem with Repulsive Manev Potential

Jagadish Singh¹, Solomon Okpanachi Omale²

¹Department of Mathematics, Faculty of Physical Sciences, Ahmadu Bello University, Zaria, Nigeria.
²Engineering and Space Systems Department, National Space Research and Development Agency (NASRDA), Obasanjo Space Centre, Abuja, Nigeria.

DOI: 10.4236/oalib.1104980                

Abstract

In this paper, we carried out a numerical study of the planar restricted four-body problem with repulsive Manev potential and perturbations in the Coriolis and centrifugal forces such that the peripherals possess Eulerian configuration. We have presented the equations of motion in the rotating frame and investigated the existence and location of the equilibrium points. We have found that there exist six equilibrium points all of which lie along the coordinate axes and shift in positions as the perturbation parameter is varied. We have also examined the linear stability of these equilibrium points and they are found unstable. The dynamical behavior of this system is also investigated using the Lyapunov Characteristic Exponents and the system is found to be chaotic.

Keywords

Repulsive-Manev Potential, Coriolis Force, Centrifugal Force, Stability, LCEs, Chaos

1. Introduction

Over the years, Mathematicians and Astronomers have been thrilled by the study of the motion of systems on n-bodies. Sir Isaac Newton pioneered the central-force and two-body problem in his work “Principia” which was first published in [1]. However, the failure of the classical gravitational law to explain the circular moon’s orbit around the earth within the frame of the inverse-square force model and other observed phenomena in the solar system dynamics such as the perihelion advances of the inner planets (for instance, mercury), got Newton to study a central-force problem given by a

potential of the type $\frac{A}{r}+\frac{B}{r^2}$. In Principia’s Book I, Article IX, Proposition

XLIV, Theorem XIV, Corollary 2, Newton showed that a central-force problem having this kind of potential leads to precessional elliptic relative orbit. That is, the trajectory of one particle considered with respect to a fixed point moves along an ellipse whose focal axis rotates in the plane of motion. Alexis Clairaut also studied this potential, but finally abandoned it in lieu of the classical potential.

There were other pre- and post-relativistic laws (such as those proposed by Hall and Newcomb) which were able to explain the phenomena of perihelion advances, but unable to justify other issues such as the secular motion of the moon’s perigee. Fortunately, the general relativity theory thrived in expounding well such phenomena in both quantitative and qualitative manner, only with the shortcomings that this powerful theory is not of much help for celestial mechanics as all attempts to formulate a meaningful relativistic n-body problem have failed to provide valuable results.

Therefore, the interest is to find a model that can maintain dynamical astronomy within the context of classical mechanics, as well as proffering justifications for the observed phenomena as offered by the relativity theory. Such a model meets the theoretical needs of celestial mechanics (by preserving the simplicity and advantages of Newtonian mechanics), and can also describe accurately the orbits coming close to collisions. By using physical principles, the Bulgarian Physicist George Manev obtained a similar model in the twenties, and proposed an alternative substitute for the relativity theory [2] [3] [4] [5] [6]. In the corresponding central force problem with unit mass for the satellite, Manev’s

potential gives $A=\mu$ and $B=\frac{3{\mu }^2}{2C^2}$, where $\mu$ is the gravitational parameter

of the two-body and C the speed of light. The Manev’s model explains the solar-system phenomena with same accuracy as relativity, but without leaving the framework of classical mechanics and it builds a bridge between the classical mechanics and the general relativity. In the recent times, researchers have taken interest in investigating the restricted few-body problem with Manev-type forces [7] [8] [9] [10] [11].

The restricted four-body problem describes the motion of an infinitesimal mass under the gravitational attraction of three massive bodies (called primaries) moving in circular orbits around their centre of mass fixed at the origin of the coordinate system. It is known that in the planar restricted four-body problem, there exist only two configurations, namely, the Eulerian (or collinear) and Lagrangian (or triangular) configurations. In the case of the later, the primaries lie at the vertices of an equilateral triangle, while in the former case, the peripherals lie on a straight line. The classical restricted four-body may be generalized to include various types of effects such as variation of the mass of the primaries, radiation pressure force, Poynting-Robertson drag, oblateness of primaries, Coriolis and centrifugal forces, etc. Several researchers such as [12] [13] [14] [15] have considered the effects of small perturbations in the Coriolis and centrifugal forces in the framework of restricted three-body problem.

In this study, our aim is to carry out a numerical investigation of the motion of an infinitesimal mass in the gravitational field of three primaries which are in Eulerian configuration under the effect of small perturbations in the Coriolis and the centrifugal forces together with the bigger primary having a repulsive Manev potential. We studied the equilibrium points, the zero velocity curves, the linear stability and the dynamical behavior of the problem with the restriction that the infinitesimal mass has no influence on the motion of the primary bodies.

2. Equations of Motion

We consider the motion of a test particle P of infinitesimal mass m under the gravitational attractions of three bodies P₁, P₂ and P₃ of masses $M_1$, $M_2$ and $M_3$ respectively, where the gravitational potential of $M_1$ is given by a Manev

potential $-\left(\frac{a}{r}+\frac{e}{r^2}\right)$ with parameter $e>0$, while the gravitational attraction

due to $M_2$ and $M_3$ is Newtonian $\left(-1/r\right)$. Also, the primaries have Euler configuration such that $M_2=M_3=\mu$ are located symmetrically with respect to the central body $M_1$, of mass $M_1=\beta \mu$, which is at the centre of masses of the system (Figure 1). In the inertial frame of reference, the peripherals $M_2$ and $M_3$ move in circular orbits about the central body $M_1$ with angular velocity $\omega$. Now, in a rotating frame Oxyz, we choose the units of the distance, mass and time such that the distance between the peripherals is unity and $G\mu =1$, where G is the gravitational constant. Let the coordinates of the infinitesimal mass and peripheral masses $M_1$, $M_2$ and $M_3$ be $\left(x,y\right)$, $\left(0,0\right)$,

$\left(\frac{1}{2},0\right)$ and $\left(-\frac{1}{2},0\right)$ respectively. As given by [16] and [17] , the peripherals

maintain their circular orbit of radius $\frac{1}{2}$ and angular velocity $\omega$ around the central body under the condition that ${\omega }^2=\Delta$, where

$\Delta =\Delta \left(\beta ,e\right)=2\left(1+4\beta -16\beta e\right)$ (1)

where $\Delta$ is a positive function, implying the parameter $e>0$ must satisfy the following sharp bound

$e<e_0=\frac{1+4\beta }{16\beta }$ (2)

here $e>0$ is admissible for a fixed value of $\beta >0$ whenever the inequality (2) is satisfied. Now, we introduce small perturbations in the Coriolis and centrifugal forces with the use of the parameters $\varphi$ and $\psi$. The unperturbed value of each is unity.

Figure 1. Sketch of the system.

The equations of motion of the infinitesimal mass under small perturbations in the Coriolis and centrifugal forces in the synodic frame can be written as:

$\begin{array}{l}\ddot{x}-2\varphi \stackrel{\dot{}}{y}={\Omega }_x\\ \ddot{y}+2\varphi \stackrel{\dot{}}{x}={\Omega }_y\end{array}\}$ (3)

where the dots denote time derivatives and the gravitational potential is given as

$\Omega =\frac{\psi \left(x^2+y^2\right)}{2}+\frac{1}{\Delta }\left[\beta \left(\frac{1}{r_1}-\frac{e}{r_1^2}\right)+\frac{1}{r_2}+\frac{1}{r_3}\right]$ (4)

with

$\begin{array}{l}{r}_1={\left(x^2+y^2\right)}^{\frac{1}{2}}\\ r_2={\left({\left(x-\frac{1}{2}\right)}^2+y^2\right)}^{\frac{1}{2}}\\ r_3={\left({\left(x+\frac{1}{2}\right)}^2+y^2\right)}^{\frac{1}{2}}\end{array}\}$ (5)

and

$\varphi =1+\epsilon ,\left|\epsilon \right|\ll 1,$

$\psi =1+{\epsilon }^{\prime },\left|{\epsilon }^{\prime }\right|\ll 1$

where $\epsilon ,{\epsilon }^{\prime }$ are small perturbations given to the Coriolis and the centrifugal forces respectively. The subscripts x and y indicate the partial derivatives of $\Omega$ with respect to x and y respectively. The system (3) possesses the energy integral

$v^2={\stackrel{\dot{}}{x}}^2+{\stackrel{\dot{}}{y}}^2=2\Omega -C$ (6)

where C is the Jacobi integral constant.

3. Location and Existence of Equilibrium Points

We investigate the existence and locations of equilibrium points of the test particle (infinitesimal mass) in this section. At these points the net force acting on the infinitesimal mass is zero. Thereby, its velocity and acceleration are both zero in the rotating frame of reference. That is, the equilibrium points satisfy $\ddot{x}=\ddot{y}=\stackrel{\dot{}}{x}=\stackrel{\dot{}}{y}=0$. It thus follows from Equation (3), that the equilibrium points are solutions of equations

$\psi x-\frac{1}{\text{Δ}}\left[\beta \left(\frac{1}{r_1^3}-\frac{2e}{r_1^4}\right)x+\frac{\left(x-\frac{1}{2}\right)}{r_2^3}+\frac{\left(x+\frac{1}{2}\right)}{r_3^3}\right]=0$ (7)

$y\left[\psi -\frac{1}{\Delta }\left[\beta \left(\frac{1}{r_1^3}-\frac{2e}{r_1^4}\right)+\frac{1}{r_2^3}+\frac{1}{r_3^3}\right]\right]=0$ (8)

we observe that Equations (7) and (8) are independent of $\varphi$. This shows that a small perturbation in the Coriolis force has no effect on the positions of equilibrium points.

Theorem (Barrabes et al., 2017); for any $\beta >0$ and admissible e, the equilibrium points of the Manev R4BP must lie on the coordinate axes.

Solving Equations (7) and (8) when the centrifugal force is unperturbed (i.e. $\psi =1$ ), we obtain six equilibrium points lying on the coordinate axes as shown in Figure 2 confirming the theorem (Barrabes et al., 2017) (Table 1).

In Figure 3, we observe that each of the equilibrium points is symmetric to another on the x and y axes respectively and the equilibrium points on the y axis form equilateral triangles with the peripherals $M_2$ and $M_3$.

In the perturbed case, we observe that as the centrifugal force perturbation parameter $\psi$ increases, the numbers of the equilibrium points does not change but the positions of the equilibrium points with respect to the peripherals change. In Figure 4, we have shown the shifting of equilibrium points.

4. Zero Velocity Surfaces

The energy integral of the problem is given by

$C=2\Omega -{\stackrel{\dot{}}{x}}^2-{\stackrel{\dot{}}{y}}^2$ (6)

where C is known as Jacobi constant. The curves of zero velocity are defined through $2\Omega =C$. This relation defines a boundary, called Hill’s surface, which separates regions where motion is allowed or forbidden. Figure 4 shows the zero velocity curves where the problem admits. The value of the Jacobi constant increases with increase in the perturbation parameter (Figure 5 and Figure 6). ( $\psi =1$, $C{L}_{1,2}=9.1949$, $\psi =1.2$, $C{L}_{1,2}=9.72996$, $\psi =1.4$, $C{L}_{1,2}=10.2055$ ), ( $\psi =1$, $C{L}_{3,4}=13.0987$, $\psi =1.2$, $C{L}_{3,4}=13.1447$, $\psi =1.4$, $C{L}_{3,4}=13.1911$ ).

5. Linear Stability of the Equilibrium Points

We examine the motion of the infinitesimal body when small displacements are given to the coordinates of the equilibrium point (x₀, y₀) under consideration. Let $\xi$ and $\eta$ be these small displacements in the coordinates such that $x=x_0+\xi$ and $y=y_0+\eta$. Then the variational equations of motion corresponding to Equations (3) are given as

Table 1. Equilibrium points with increase in centrifugal force.

Figure 2. Showing the six equilibrium points each lying on the coordinate axes when $\beta =10$ and $e=0.25$.

Figure 3. Symmetric equilibrium points.

Figure 4. Shifts in position of equilibrium points with increase in the centrifugal force.

Figure 5. The zero velocity curves when ( $\psi =1$, $C{L}_{1,2}=9.1949$, $\psi =1.2$, $C{L}_{1,2}=9.72996$, $\psi =1.4$, $C{L}_{1,2}=10.2055$ ). With the increase in the energy constant the infinitesimal mass is trapped within the region of each of the primaries.

Figure 6. Zero velocity curves when ( $\psi =1$, $C{L}_{3,4}=13.0987$, $\psi =1.2$, $C{L}_{3,4}=13.1447$, $\psi =1.4$, $C{L}_{3,4}=13.1911$ ). The white region is the region of permissible motion for the infinitesimal fourth body.

$\begin{array}{l}\ddot{\xi }-2\varnothing \stackrel{\dot{}}{\eta }=\xi \left({\Omega }_{xx}^0\right)+\eta \left({\Omega }_{xy}^0\right)\\ \ddot{\xi }-2\varnothing \stackrel{\dot{}}{\eta }=\xi \left({\Omega }_{xx}^0\right)+\eta \left({\Omega }_{xy}^0\right)\end{array}\}$ (9)

where the superscript 0 indicates that the values are evaluated at the equilibrium point (x₀, y₀), the subscripts represent the second partial derivatives and the dots signify the derivatives with respect to the actual time t. Here the linear terms in $\xi$ and $\eta$ are only considered.

Let the trial solutions of Equations (9) be

$\xi =P{\text{e}}^{\lambda t},\eta =Q{\text{e}}^{\lambda t}$ (10)

where $P,Q$ are constants and $\lambda$ is a parameter. Then the characteristic equation of the system (10) can be written as

${\lambda }^4+a{\lambda }^2+b=0$ (11)

With

$a=4{\varphi }^2-{\Omega }_{xx}^0-{\Omega }_{yy}^0$

$b={\Omega }_{xx}^0{\Omega }_{yy}^0-{\left({\Omega }_{xy}^0\right)}^2$

$\begin{array}{l}{\Omega }_{xx}^0=\psi +\frac{1}{\Delta }[\beta \left(-\frac{1}{r_{10}^3}+\frac{3x^2}{r_{10}^5}+\frac{2e}{r_{10}^4}-\frac{8e{x}^2}{r_{10}^6}\right)+\frac{3{\left(x-1/2\right)}^2}{r_{20}^5}\\ -\frac{1}{r_{20}^3}+\frac{3{\left(x+1/2\right)}^2}{r_{30}^5}-\frac{1}{r_{30}^3}]\end{array}$ (12)

${\Omega }_{yy}^0=\psi +\frac{1}{\Delta }\left[\beta \left(-\frac{1}{r_{10}^3}+\frac{3y^2}{r_{10}^5}+\frac{2e}{r_{10}^4}-\frac{8e{y}^2}{r_{10}^6}\right)+\frac{3y^2}{r_{20}^5}-\frac{1}{r_{20}^3}+\frac{3y^2}{r_{30}^5}-\frac{1}{r_{30}^3}\right]$ (13)

${\Omega }_{xy}^0=\frac{1}{\Delta }\left[\beta \left(\frac{3xy}{r_{10}^5}-\frac{8exy}{r_{10}^6}\right)+\frac{3\left(x-\frac{1}{2}\right)y}{r_{20}^5}+\frac{3\left(x+\frac{1}{2}\right)y}{r_{30}^5}\right]$ (14)

$\begin{array}{l}{r}_{10}={\left(x_0^2+y_0^2\right)}^{\frac{1}{2}}\\ r_{20}={\left({\left(x_0-\frac{1}{2}\right)}^2+y_0^2\right)}^{\frac{1}{2}}\\ r_{30}={\left({\left(x_0+\frac{1}{2}\right)}^2+y_0^2\right)}^{\frac{1}{2}}\end{array}$

The four roots of the characteristic Equation (11) play an important role in the determination of stability of the equilibrium points. An equilibrium point under consideration will be stable if the Equation (11) has all four purely imaginary roots or has four complex roots with each of them having negative real part. This is equivalent to saying that the following system of inequalities must be simultaneously satisfied (Table 2).

Table 2. Stability of Equilibrium points.

${\left(4{\varphi }^2-{\Omega }_{xx}^0-{\Omega }_{yy}^0\right)}^2-4\left({\Omega }_{xx}^0{\Omega }_{yy}^0-{\left({\Omega }_{xy}^0\right)}^2\right)>0$

$\left(4{\varphi }^2-{\Omega }_{xx}^0-{\Omega }_{yy}^0\right)>0$

$\left({\Omega }_{xx}^0{\Omega }_{yy}^0-{\left({\Omega }_{xy}^0\right)}^2\right)>0$

We have computed the characteristic roots of Equation (11) as perturbation parameters $\varphi$ and $\psi$ increase and we observe that the equilibrium points L_i (i = 1, 2, 3, 4, 5, 6) are unstable.

6. Dynamic Behaviour of the System

The Lyapunov Characteristic Exponents (LCEs) measure the average rate of convergence or divergence of orbits starting from nearby orbits. It tells whether or not two points in the phase space of a dynamical system that are initially very close will remain close as the motion of the system proceed. It is used as a tool to describe the behaviour of the dynamical systems. It is employed to determine the existence of chaos or regularity of the orbits (for example see [18] ). In the applicable sense, the exponential divergence of the orbits connotes impossibility to predict the system, so according to [19] any system with at least one positive Lyapunov exponent is chaotic. Therefore, LCEs can be used to analyze the stability of limit sets and to check sensitive dependence on initial conditions, that is, the presence of chaotic attractors (Figure 7).

We have computed numerically the first order LCEs and plotted the graphs (LCEs vs Steps) with the help of Mathematica package developed by Sandri [20]. It is a customary practice to refer to the Maximal Lyapunov Exponent (MLE), because it determines a notion of predictability for a dynamical system. We find that the system is chaotic because the LCEs [0.965343, 0.965343, −0.965343, −0.965343] contain two positive exponents.

7. Discussion and Conclusion

We have studied the existence, location, stability and dynamical behavior of the equilibrium points of an infinitesimal mass under small perturbations in the Coriolis and centrifugal forces in the restricted four-body problem when the peripherals have Eulerian configuration with a repulsive Manev potential. We have expressed in the rotating coordinate system the equations governing the motion of the infinitesimal mass, and using the energy integral we have determined the region of permissible motion by the zero velocity curves. We observe that the Coriolis force has no effect on both the location of the equilibrium points and the zero velocity curves, but a small perturbation in the centrifugal force affects both the positions of the equilibrium points and the zero velocity curves. We have found six equilibrium points all located on the coordinate axes that verify numerically the theorem of Barrabes et al. (2017). We also observe that in addition to perturbations in the Coriolis and centrifugal forces which cause the orbit of the infinitesimal body to shrink, the equilibrium points L_i (i = 1, 2, 3, 4, 5, 6)

Figure 7. The Lyapunov Characteristic Exponents of the system.

are unstable. With the aid of Mathematica package, we also computed the LCEs of the system and found that the system is chaotic, because two of its exponents are positive.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Newton, I. (1999) The Principia: Mathematical Principles of Natural Philosophy. Uni-versity of California Press, ?Oakland, CA.
[2]	Maneff, G. (1924) La gravitation et le principe de l’egalie de l’action et de la reaction. Comptes Rendus de l’Académie des Sciences de Paris, 178, 2159-2161.
[3]	Maneff, G. (1925) Die Gravitation und das Prinzip von Wirkung und Gegenwirkung. Zeitschrift für Physik, 31, 786-802. https://doi.org/10.1007/BF02980633
[4]	Maneff, G. (1929) Die Masse der Feldenergie und die Gravitation. Astronomische Nachrichten, 236, 401-406. https://doi.org/10.1002/asna.19292362402
[5]	Maneff, G. (1930) La gravitation et l’energie au zero. Comptes Rendus de l’Académie des Sciences de Paris, 190, 1374-1377.
[6]	Blaga, C. (2015) Prescessing Orbits, Central Forces and Manev Po-tential. In: Gerdjikov, V. and Tsetkov, M., Eds., Prof. G. Manev’s Legacy in Contempo-rary Aspects of Astronomy, Theoretical and Gravitational Physics, Heron Press Ltd., Chicago, IL, 134-139.
[7]	Ivanov, R. and Prodanov, E. (2005) Manev Potential and General Relativity. In: Gerdjikov, V. and Tsetkov, M., Eds., Prof. G. Manev’s Legacy in Contemporary Aspects of Astronomy, Theoretical and Gravitational Physics, Heron Press Ltd., Chicago, IL, 148-154.
[8]	Haranas, I. and Mioc, V. (2009) Manev Potential and Satellite Orbits. Romanian Astronomical Journal, 19, 153-166.
[9]	Kirk, S., Haranas, I. and Gkigkitzis, I. (2013) Satellite Motion in a Manev Potential with Drag. Astrophysics and Space Science, 344, 313-320. https://doi.org/10.1007/s10509-012-1330-0
[10]	Blaga, C. (2015) Stability in Sense of Lyapunov of Circular Orbits in Manev Potential. Romanian Astronomical Journal, 25, 233-240.
[11]	Barrabes, E., Cors, J. and Vidal, C. (2017) Spatial Collinear Restricted Four Body Problem with Repulsive Manev Potential. Celestial Mechanics and Dynamical Astronomy, 129, 153-176. https://doi.org/10.1007/s10569-017-9771-y
[12]	Bhatnager, K.P. and Hallan, P.P. (1978) Effect of Perturbations in Coriolis and Centrifugal Forces on the Stability of Liberation Points in the Restricted Problem. Celestial mechanics, 18, 105-112. https://doi.org/10.1007/BF01228710
[13]	Singh, J. and Vincent, A.E. (2015) Effect of Perturbations in the Coriolis and Centrifugal Forces on the Stability of Equilibrium Points in the Restricted Four-Body Problem. Few-Body Systems, 56, 713-723. https://doi.org/10.1007/s00601-015-1019-3
[14]	Abdul Raheem, A. and Singh, J. (2006) Combined Effects of Perturbations, Radiation, and Oblateness on Stability of Equilibrium Points in the Restricted Three-Body Problem. The Astronomical Journal, 131, 1880-1885. https://doi.org/10.1086/499300
[15]	Abouelmagad, I.A., Asiri, H.M. and Sharaf, M.A. (2013) The Effect of Oblateness in the Perturbed Restricted Three-Body Problem. Meccanica, 48, 2479-2490. https://doi.org/10.1007/s11012-013-9762-3
[16]	Arribas, M., Elipe, A. and Kalvouridis, T. (2007) Periodic Solutions in the Planar (n+1) Ring Problem with Ob-lateness. Journal of Guidance, and Dynamics, 30, 1640-1648. https://doi.org/10.2514/1.29524
[17]	Fakis, D.G. and Kalvouridis, T.J. (2013) Dy-namics of a Small Body under the Action of a Maxwell Ring-Type N-Body System with a Spheroidal Central Body. Celestial Mechanics and Dynamical Astronomy, 116, 229-240. https://doi.org/10.1007/s10569-013-9484-9
[18]	Kumari, R. and Kushvah, B.S. (2013) Equilibrium Points and Zero Velocity Surfaces in the Restricted Four-Body Problem with Solar Wind Drag. Astrophysics and Space Science, 344, 347-359. https://doi.org/10.1007/s10509-012-1340-y
[19]	Dubeibe, F.L. and Bermudez-Almanza, L.D. (2013) Optimal Conditions for the Numerical Calculation of the Largest Lyapunov Exponent for Systems of Ordinary Differential Equations. International Journal of Modern Physics C, 25, Article ID: 1450024. https://doi.org/10.1142/S0129183114500247
[20]	Sandri, M. (1995) Numerical Calculation of Lyapunov Exponents. University of Verona, Italy. http://www.msandri.it/docs/lce.m