universe Article Approximate Analytical Periodic Solutions to the Restricted Three-Body Problem with Perturbation, Oblateness, Radiation and Varying Mass Fabao Gao 1,2, * and Yongqing Wang 1 1 School of Mathematical Science, Yangzhou University, Yangzhou 225002, China; yqwangyzu@gmail.com 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain * Correspondence: fbgao@yzu.edu.cn  Received: 17 June 2020; Accepted: 2 August 2020; Published: 4 August 2020  Abstract: Against the background of a restricted three-body problem consisting of a supergiant eclipsing binary system, the two primaries are composed of a pair of bright oblate stars whose mass changes with time. The zero-velocity surface and curve of the problem are numerically studied to describe the third body’s motion area, and the corresponding five libration points are obtained. Moreover, the effect of small perturbations, Coriolis and centrifugal forces, radiative pressure, and the oblateness and mass parameters of the two primaries on the third body’s dynamic behavior is discussed through the bifurcation diagram. Furthermore, the second- and third-order approximate analytical periodic solutions around the collinear solution point L3 in two-dimensional plane and three-dimensional spaces are presented by using the Lindstedt-Poincaré perturbation method. Keywords: periodic solution; bifurcation analysis; three-body problem; Lindstedt-Poincaré perturbation 1. Introduction It is well known that Newton’s famous proposition “Principle” took a decisive step in our understanding of the universe. The origin of the universe is inseparable from the three-body dynamics, which has always existed since the universe’s birth. It can be said that there are no three-body dynamics; it is impossible to create a universe [1]. The three-body interactions in the early universe help people understand the black hole MACHO (abbreviation of MAssive Compact Halo Object) binary formation and the detectability of gravitational waves in black hole MACHO binaries coalesce [2,3]. In the theoretical framework of scalar-tensor gravitation, Zhou et al. [4] found that there is a collinear solution to the three-body problem in the presence of a scalar field, and studied the effect of the scalar field on this solution and the positions of Lagrange points through numerical examples. For the three-body problem in the spherical universe, their perturbation theory analysis showed that the rate of precession of two small and nearly circular solutions of identical particles is proportional to the square root of their initial distance and inversely proportional to the square of the radius of the universe [5]. In addition, the three-body problem is also widely used in the evolution of binary systems [6] and the dynamic analysis of binary asteroids [7,8], as well as other fields in the universe such as dark matter, galaxies, GW170817 (GW is short for gravitational wave), and Mukhanov-Sasaki Hamiltonian dynamics, and so forth (see References [9–13] for more information). The mathematician Poincaré believed periodic solutions to be the unique avenue for solving the three-body problem. He stated that “what makes these solutions so precious to us, is that they are, so to say, the only opening through which we can try to enter a place.” There is too much literature in this area, mainly involving three aspects—qualitative analysis, approximate quantitative solution, and computer-aided numerical simulation (see References [14–32], and the references therein). Universe 2020, 6, 110; doi:10.3390/universe6080110 www.mdpi.com/journal/universe Universe 2020, 6, 110 2 of 25 When considering multiple perturbation factors, Singh et al. did a great deal of work on the libration points and their stability of various restricted three-body models (see References [33–41]). These factors mainly include the oblateness of the primaries [35–38,40,41], the oblateness of the third body [37], and the small perturbation in centrifugal force and Coriolis force [35–41], the radiation pressure of the primaries [33–41], the changing mass of the primaries over time based on the Meshcherskii law [34–36], the changing mass of the third body governed by Jeans law [42], the triaxiality of the primaries [39]. Basic references on oblateness, Coriolis and centrifugal forces, and radiation pressure can be found in References [43–47]. In addition, some other researchers also discussed the periodic solutions and libration points of the restricted three-body problem (R3BP) under specific perturbations. For example, Abdullah et al. [48] investigated the locations and stability of libration points in the circular R3BP with Coriolis and centrifugal forces under the assumption that the two primaries and the infinitesimal body with variable masses. Under the assumption that the more massive primary has triaxiality, Abouelmagd et al. [49] obtained the approximate periodic solution near the collinear libration point of R3BP to the second-order. Since the star’s isotropic radiation or absorption is one of the factors that cause the mass of the celestial bodies to change with time, Bekov [50] considered the R3BP of variable masses, in which the motion of the two primaries was determined by the Gylden-Meshcherskii problem and found the collinear solutions, as well as the spatial solution and its domains of existence. When the more extensive primary is an oblate sphere, Zotos [51] numerically studied the dynamic behavior of the circular R3BP under the initial conditions of the solution are bounded, escaping and collisional, and found that the oblateness factor has a significant influence on the characteristics of the solutions. Based on the model of the binary system proposed in Reference [39], Gao and Wang [6] continued to study the analytical approximate periodic solutions around the collinear libration points. They examined the influence of the small perturbation in Coriolis and centrifugal forces, the triaxiality, and radiation pressure of the primaries on the third body’s dynamic behavior through the bifurcation diagram. In this paper, the effects of small perturbations, Coriolis and centrifugal forces, radiation pressure, and the oblateness and variable mass of the primaries on the third body’s dynamic behavior will be discussed numerically through bifurcation diagram in the context of an R3BP consisting of a supergiant eclipsing binary system. Here we take into account the oblateness, Coriolis and centrifugal forces together because the oblateness of the primaries will cause both Coriolis force and centrifugal force to increase (see References [52,53]). By Lindstedt-Poincaré perturbation method, the approximate analytical solutions of the perturbed third body in two and three dimensions near the collinear libration points will be obtained. This paper is organized as follows—the motion equations in an autonomous system will be given in Section 2. The geometry structure of the zero velocity surfaces and curves will be illustrated in Section 3. Moreover, the bifurcation diagrams of the different perturbation parameters will be shown in Section 4. Furthermore, the periodic planar and spatial solutions near the collinear libration points will be obtained by means of Lindstedt-Poincaré in Section 5. In the last Section 6, conclusions will be stated. 2. Equations of Motion In a barycentric coordinate Oxyz, consider that two oblate variable-mass radiating primaries move at an angular velocity ω in a circular orbit around their center of mass. The plane of motion of the primaries (located on the x-axis) coincides with the x-y plane. The third body of infinitesimal mass that moves under the gravity of the primaries but does not affect their motion, then the dynamic governing equations (see Reference [35] or Reference [36] for more details) can be written as follows Universe 2020, 6, 110 3 of 25 d2 x dy dω µ ( t ) q1 ( x − x1 ) µ2 ( t ) q2 ( x − x2 ) − 2ω α − αy = ω 2 xβ − 1 − dt2 dt dt r13 r23 3µ1 (t)q1 A1 (t) ( x − x1 ) 3µ2 (t)q2 A2 (t) ( x − x2 ) − − , 2r15 2r25 (1) d2 y dx dω µ (t)q y µ2 (t)q2 y 3µ1 (t)q1 A1 (t)y 3µ2 (t)q2 A2 (t)y + 2ω α + αx = ω 2 yβ − 1 3 1 − − − , dt2 dt dt r1 r23 2r15 2r25 d2 z µ (t)q z µ (t)q z 3µ (t)q1 A1 (t)z 3µ2 (t)q2 A2 (t)z =− 1 3 1 − 2 3 2 − 1 − , dt2 r1 r2 2r15 2r25 where ri2 = ( x − xi )2 + y2 + z2 (i = 1, 2) denotes the distance between the ith primary and the third body, and the two primaries are placed on ( x1 , 0, 0) and ( x2 , 0, 0), respectively. For the distance r between the primaries, x1 = −µ2 (t)r/µ(t), x2 = µ1 (t)r/µ(t), µ(t) = µ1 (t) + µ2 (t), µi (t) = f mi (t), mi (t) represents the mass of the ith primary, f is the gravitational constant. In addition, the parameters qi and Ai (t) stand for the radiation factor and the varying oblateness of the ith primary, respectively. The Coriolis force α and centrifugal force β satisfy α = 1 + ε, β = 1 + ε0 , here |ε| and |ε0 | are the corresponding small perturbations with values far less than 1. Furthermore, it can be followed from the reference [36] that the angular velocity holds the following expression   23µ(t) 2 A1 ( t ) + A2 ( t ) ω = + , (2) 2κr3 (t) 3 r2 (t) where κ is an arbitrary dimensionless constant of a particular integral rµ(t) = κC2 in Gyldén-Meshcherskii problem (see References [54–56] for more information), it represents an arbitrary sum of the masses of the primaries, C is a constant of the area integral. Note that if the masses of the primaries are constant and the sum κ = 1, then µ(t) = 1 by unitizing the gravitational constant. In addition, if the distance r between the primaries is 1, and without considering the oblateness of the second primary, that is, A2 (t) = 0. Then Equation (2) is simplified to ω 2 = 1 + 3A1 (t)/2, which is the same as the mean motion expression in Reference [47]. We now transform Equations (1) into the autonomous system (u, v, w, τ ) by the Meshcherskii’s transformation (see Reference [35] or Reference [57]) x = u(τ ) R(t), y = v(τ ) R(t), z = w(τ ) R(t), dt/dτ = R2 (t), ri = di R(t), (3) p where di is the distance between the third body and the ith primary, R(t) = φt2 + 2ϕt + γ, φ, ϕ and γ are constants, the unified Meshcherskii law µ(t) = µ0 /R(t), µ1 (t) = µ10 /R(t), µ2 (t) = µ20 /R(t), (4) where µ0 , µ10 , µ20 are constants, and the particular solutions of Gyldén-Meshcherskii problem, x1 = u1 R(t), x2 = u2 R(t), C = d212 ω0 , r = d12 R(t), ω = ω0 /R2 (t), (5) where d12 is the distance between the two primaries at initial time. Applying transformation Ai (t) =   αi R2 (t), here αi = d2Ei − d2Pi 5d212 (i = 1, 2) refers to the oblateness of the ith primary at initial  time, d Ei and d Pi are the equatorial radii and the polar radii of the ith primary, respectively. Similarly, following Singh and Leke [35], unify the mass and distance at the initial time and introduce the mass parameter µm which satisfies µ10 /µ0 = 1 − µm , µ20 /µ0 = µm , µm ∈ (0, 1/2] stands for the ratio of the smaller mass to the total mass of the two primaries, then µ0 = κ, d12 = 1. Therefore, we obtain " # 2 α1 R2 ( t ) + α2 R2 ( t )   2 3µ0 1 3 ω = + = 1 + ( α + α 2 ) . (6) 2κd312 R4 (t) 3 d212 R2 (t) R4 ( t ) 2 1 Universe 2020, 6, 110 4 of 25 It follows from Equation (5) that ω02 = 1 + 3(α1 + α2 )/2, αi  1 (i = 1, 2). Then Equations (1) are reduced to the following form ∂Ω ü − 2αω0 v̇ = , ∂u ∂Ω v̈ + 2αω0 u̇ = , (7) ∂v ∂Ω ẅ = , ∂w with the potential function " # ω02 ( β + κ − 1)  2  ω 2 (κ − 1) w2 q (1 − µ m ) q µm q α (1 − µ m ) q α µm Ω= u + v2 + 0 +κ 1 + 2 + 1 1 3 + 2 23 , (8) 2 2 d1 d2 2d1 2d2 q q 2 where d1 = (u + µm ) + v2 + w2 and d2 = (u + µm − 1)2 + v2 + w2 are the distances from the third body to the two primaries, respectively. 3. Zero Velocity Surfaces and Curves The level surface of Equation (7) is an energy surface. The projection of this surface onto position space is named Hill’s region, and its boundary is the zero-velocity surface, on which the velocity of the third body is zero, so the third body can only move inside the zero-velocity surface. The intersection of this surface with the x-y plane is known as the zero-velocity curve. Note that Equation (7) admits a first integral   2Ω − u̇2 + v̇2 + ẇ2 = C, (9) where C is the Jacobi constant and u̇2 + v̇2 + ẇ2 denotes the velocity of the relative motion. The geometric structure of 2Ω = C in the two-dimensional plane and three-dimensional space is shown as a series of curves and surfaces for the different Jacobi constants C. The zero-velocity surface (or curve) with perturbations helps us to analyze the third body’s dynamic behavior under the influence of a specific perturbation factor, which is different from that of classical R3BP. For example, when q1 = q2 = 1, α1 = α2 = 0 and α = β = 1, the zero-velocity surface is different from the classic R3BP due to the presence of κ. For the parameter values presented in Table 1, Figure 1 shows the four zero-velocity surfaces when κ takes 0.5, 1, 2 and 4, respectively. It is easy to find that the structure of the zero-velocity surface of Equation (9) is more open as the value of κ becomes smaller and smaller. Because the geometry of the zero-velocity surface (or curve) will vary with the value of the Jacobi constant C and κ, we can numerically analyze the zero-velocity surface (or curve) and then find the libration points of Equation (7). Without loss of generality, now we only study the surfaces and curves for the case κ = 1. For 2.82017 ≤ C ≤ 4.16993, the forbidden region of the third body changes with the Jacobi constant C as shown in Figures 2 and 3. With the decrease of C, five regular libration points L1 , L2 , · · · , L5 are discovered gradually and the forbidden region of motion of the third body is shrinking. More precisely, the third body can only fly in the two ellipsoidal areas around the two primaries when C = 4.16993, and the two ellipsoidal motion areas are connected together through the “tangent point” (marked as librarian point L1 ), and the third body can pass from one primary to the other through L1 when C < 4.16993, but it cannot fly out of these two ellipsoidal areas. When C = 3.59308, the peanut-shaped flight area is connected to the outer space through a new “tangent point” L2 , and when the value of C continues to decrease, the third body can only enter the outer space in this direction. When C drops to 3.57940, we will discover another new “tangent point” in the opposite direction, which is the libration point L3 . At this time the third body can also fly from this direction to outer space when the value of C continues to decrease. For the Jacobi constant C = 2.82017, Universe 2020, 6, 110 5 of 25 the symmetrical libration points L4 and L5 appear at the same time. Now the third object can escape into outer space in any direction around the two primaries. Therefore, as a summary of Figures 2 and 3, the corresponding five planar librations points can be identified as follows L1 (0.018995, 0), L2 (1.1989, 0), L3 (−1.19261, 0), L4 (0.01424987, 0.851953), L5 (0.01424987, −0.851953). Besides, there are several methods for determining the libration points. For example, the total number of collinear libration points for this problem can also be determined based on the theory of topological degree [58]. Figure 1. Zero velocity surfaces for κ = 0.5, 1, 2 and 4, respectively. (a) (b) (c) (d) Figure 2. Zero-velocity surfaces of the third body for κ = 1: (a) C = 4.16993, (b) C = 3.59308, (c) C = 3.57940, (d) C = 2.82017. Universe 2020, 6, 110 6 of 25 Table 1. Physical parameters in Equation (7). Parameters µm q1 q2 α1 α2 α β Values 0.48785 0.9988 0.9985 0.024 0.02 1.001 1.002 2 2 0 0 -2 -2 -2 0 2 -2 0 2 (a) (b) 2 2 0 0 -2 -2 -2 0 2 -2 0 2 (c) (d) Figure 3. Zero-velocity curves of the third body: (a) C = 4.16993, (b) C = 3.59308, (c) C = 3.57940, (d) C = 2.82017. 4. Bifurcation Analysis We note that the three-body problem has a high degree of complexity due to its dynamic equations’ extremely strong nonlinearity. Any subtle changes in different parameters may substantially change the third body’s dynamic behavior. Thus, we will analyze the effect of different perturbation parameters on the third body’s dynamic behavior through the bifurcation diagram. After a large number of numerical simulations, we select u0 = (−0.08, 0.001, 0.2, 0.002, −1.6, 0.001) as the initial value of the iteration. Combining the conditions in section 2, we restrict the Coriolis force α and centrifugal force β to [1, 1.1] in the α-β-u frame, Figure 4a,b show the influence of the perturbation of α and β on the variables u-v, respectively, and reflect the corresponding chaotic characteristics. However, the joint perturbation of α and β shows that the amplitude of the third body in u direction no longer changes when 1.055 < α < 1.1 (see Figure 4c). Universe 2020, 6, 110 7 of 25 (a) (b) (c) Figure 4. (a–c) Bifurcation diagrams of Coriolis force α and centrifugal force β in the frames αuv, βuv and αβu, respectively, when κ = 1, µm = 0.48785, q1 = 0.9988, q2 = 0.9985, α1 = 0.024, α2 = 0.02. Universe 2020, 6, 110 8 of 25 Consider the perturbation range of the oblateness coefficients of the two primaries is [0, 0.1]. As shown in Figure 5a,b, as the oblateness coefficient α1 of the first primary increases, the effect of this perturbation factor on the third body’s amplitude in u direction gradually increases, but the effect on the amplitude in v direction gradually decreases. When the oblateness coefficient α2 of the second primary changes, the third body exhibits apparent chaotic behavior. However, the third body keeps moving periodically in u direction under the joint action of α1 and α2 (see Figure 5c). This indicates that the oblateness coefficients play a certain degree of “correction” to the two primaries’ dynamic behavior so that they can maintain regular motion as much as possible. For the parameter κ, which represents the arbitrary sum of the primaries’ masses, we only show the bifurcation diagrams whose value is between 0 and 1. Figure 6a–c clearly show that no matter how the value of κ changes, it has a considerable influence on the third body’s dynamic behavior. The third body reflects the intricate, chaotic properties in the frames κ-u-v, κ-u-w, and κ-v-w. It is consistent with our common sense because the masses of the primaries are closely related to their gravitation, and this plays a significant role in determining the behavior of the third body in the sphere of gravity. Next, we show the bifurcation diagrams of the mass parameter µm in the interval (0, 0.5] in Figure 7. It is easy to find that when µm ’s value is small, the mass of one primary is enormous compared to the other. It has a significant impact on the third body’s amplitude in the three directions of u, v, and w. When µm = 0.07, the third body’s amplitude in u direction seems to be stable from its projection (see Figure 7a,b), but it is still apparent in the three-dimensional space. However, when the value of µm gradually increases, that is, the difference between the masses of the two primaries is getting smaller and smaller, its influence on the third body’s amplitude in the v direction is gradually more considerable than that in u and w directions (see Figure 7b,c). Similarly, we examine the effect of the radiation factors q1 and q2 on the third body’s dynamic behavior at the interval [0, 1]. Figure 8a,b demonstrate the influence of the radiation factors of the first and second primaries, respectively, and Figure 8c reflects the impact of the two radiation factors acting simultaneously. As shown in Figure 8a,b, the radiation factor q1 or q2 has almost the same influence on the third body’s movement in u and v directions. However, the influence of q1 on the third body’s amplitude is greater than that of q2 , especially when q1 < 0.57 or 0.67 < q1 < 0.85. In contrast, the impact of q2 is much more moderate. Figure 8c shows that when q1 ≤ 0.0202, regardless of how q2 changes, the amplitude of the third body changes little in the u direction. In other words, the two radiations from the two primaries in the current state have almost little effect on the dynamic behavior of the third body in the u direction, which may be due to the weak radiation of the first primary. However, when q1 ∈ [0.0202, 0.3434], the third body’s amplitude has a significant change, even if q2 is very small. When the value of q1 is greater than 0.3434, the third body exhibits approximately symmetrical motion in a larger area under the joint action of q1 and q2 . Universe 2020, 6, 110 9 of 25 0.1 0 −0.1 −0.2 v −0.3 −0.4 −0.5 0 −0.2 0.1 0.08 −0.4 0.06 −0.6 0.04 0.02 u −0.8 0 α1 (a) (b) 0.01 0.005 0 u −0.005 −0.01 0.1 0.08 0.1 0.06 0.08 0.04 0.06 0.04 0.02 0.02 α2 0 0 α1 (c) Figure 5. (a–c) Bifurcation diagrams of oblateness coefficient of the ith (i = 1, 2) primary in the frames α1 uv, α2 uv and α1 α2 u, respectively, when κ = 1, µm = 0.48785, q1 = 0.9988, q2 = 0.9985, α = 1.001, β = 1.002. Universe 2020, 6, 110 10 of 25 1.5 1 0.5 v 0 −0.5 −1 1 0.5 1 0 0.8 −0.5 0.6 −1 0.4 −1.5 0.2 u −2 0 κ (a) 0.5 0 w −0.5 0.5 0 1 −0.5 0.8 −1 0.6 0.4 −1.5 0.2 u −2 0 κ (b) 0.6 0.4 0.2 w 0 −0.2 −0.4 1.5 1 1 0.5 0.8 0 0.6 0.4 −0.5 0.2 v −1 0 κ (c) Figure 6. (a–c) Bifurcation diagrams of the arbitrary sum κ of the masses of the primaries in the frames κuv, κuw and κvw, respectively, when µm = 0.48785, q1 = 0.9988, q2 = 0.9985, α1 = 0.024, α2 = 0.02, α = 1.001, β = 1.002. Universe 2020, 6, 110 11 of 25 0.05 0 −0.05 v −0.1 0 0.5 −0.05 0.4 0.3 −0.1 0.2 0.1 u 0 µm (a) 0.25 0.2 0.15 w 0.1 0.05 0 0 0.5 −0.05 0.4 0.3 −0.1 0.2 0.1 u 0 µm (b) 0.25 0.2 0.15 w 0.1 0.05 0 0.05 0 0.5 0.4 −0.05 0.3 −0.1 0.2 0.1 v 0 µm (c) Figure 7. (a–c) Bifurcation diagrams of the mass parameter µm in the frames µm uv, µm uw and µm vw, respectively, when κ = 1, q1 = 0.9988, q2 = 0.9985, α1 = 0.024, α2 = 0.02, α = 1.001, β = 1.002. Universe 2020, 6, 110 12 of 25 15 10 5 v 0 −5 −10 10 5 1 0 0.8 −5 0.6 0.4 −10 0.2 u −15 0 q1 (a) 1.5 1 0.5 v 0 −0.5 −1 1.5 1 1 0.5 0.8 0 0.6 −0.5 0.4 −1 0.2 u −1.5 0 q2 (b) (c) Figure 8. (a–c) Bifurcation diagrams of radiation factor of the ith (i = 1, 2) primary in the frames q1 uv, q2 uv and q1 q2 u, respectively, when κ = 1, µm = 0.48785, α1 = 0.024, α2 = 0.02, α = 1.001, β = 1.002. Universe 2020, 6, 110 13 of 25 5. Periodic Solutions near the Collinear Libration Points 5.1. Expansion of Two-Dimensional Dynamic Equations The third body’s equations of motion in the u-v plane are ü − 2αω0 v̇ = Uu , (10) v̈ + 2αω0 u̇ = Uv , where the potential function " # ω02 ( β + κ − 1)  2 2  q1 (1 − µ m ) q2 µ m q1 α1 (1 − µ m ) q2 α2 µ m U= u +v +κ + + + , (11) 2 da db 2d3a 2d3b q q where d a = (u + µm )2 + v2 and db = (u + µm − 1)2 + v2 are the distances from the third body to two primaries respectively. Then the right hand side of Equation (10) can be written as q1 κ (1 − µ m ) ( u + µ m ) q κµm (u + µm − 1) 3q κα (1 − µm ) (u + µm ) Uu = ω02 ( β + κ − 1) u − 3 − 2 3 − 1 1 da db 2d5a 3q2 κα2 µm (u + µm − 1) − , (12) 2d5b q1 κ (1 − µ m ) v q κµm v 3q κα (1 − µm ) v 3q κα µm v Uv = ω02 ( β + κ − 1) v − − 2 3 − 1 1 5 − 2 25 . d3a db 2d a 2db Now we study the periodic solutions of the third body around the collinear libration points Li (i = 1, 2, 3), and u Li is the specific coordinate corresponding to Li . By letting u = u Li + ξ and v = η, then Equation (10) are reduced to the following form ξ̈ − 2αω0 η̇ = Uξ , (13) η̈ + 2αω0 ξ̇ = Uη . Expanding the right hand side of Equation (13) to the second-order yields ξ̈ − 2αω0 η̇ = K1 ξ + K2 ξ 2 + K3 η 2 , (14) η̈ + 2αω0 ξ̇ = P1 η + P2 ξη, and to the third-order in the same way, then we obtain ξ̈ − 2αω0 η̇ = K1 ξ + K2 ξ 2 + K3 η 2 + K4 ξ 3 + K5 ξη 2 , (15) η̈ + 2αω0 ξ̇ = P1 η + P2 ξη + P3 η 3 + P4 ξ 2 η, where corresponding coefficients can be found in Appendix A. Two-Dimensional Periodic Solutions Assuming that the periodic solutions for Equation (14) with respect to orbital parameter e (|e|  1) have the following form ξ = ξ 1 e + ξ 2 e2 , (16) η = η1 e + η2 e 2 . Substituting Equation (16) into Equation (14) and comparing the terms with the same power of e, we have Universe 2020, 6, 110 14 of 25 ξ̈ 1 − 2αω0 η̇1 − K1 ξ 1 = 0, (17) η̈1 + 2αω0 ξ̇ 1 − P1 η1 = 0, and ξ̈ 2 − 2αω0 η̇2 − K3 η12 − K1 ξ 2 − K2 ξ 12 = 0, (18) η̈2 + 2αω0 ξ̇ 2 − P1 η2 − P2 ξ 1 η1 = 0. By using the Lindstedt-Poincaré perturbation method (see Chapter 5 in [59] for more details), the periodic solution of Equation (17) can be written as ξ 1 = a1 cos (vt) + b1 sin (vt) , (19) η1 = a2 cos (vt) + b2 sin (vt) , with period T = 2π/v. Substituting Equation (19) into Equation (17) and comparing the coefficients of sin (vt) and cos (vt), then we get   K1 + v 2 a1 + 2αω0 vb2 = 0,   (20) 2αω0 va1 + P1 + v 2 b2 = 0, and   K1 + v 2 b1 − 2αω0 va2 = 0,   (21) − 2αω0 vb1 + P1 + v 2 a2 = 0. According to the homogeneous linear Equations (20) and (21), in order to get a nonzero solution, the corresponding determinant satisfies   v 4 + K1 + P1 − 4α2 ω02 v 2 + K1 P1 = 0, (22) then we set a1 = 1, b1 = 0, and substitute them into Equations (20) and (21) respectively, such that a2 = 0 and b2 = −2αω0 v P1 + v 2 , so the periodic solution of Equation (17) is   ξ 1 = cos (vt) , (23) η1 = b2 sin (vt) . Similarly, we write the periodic solution of Equation (18) as ξ 2 = a3 + a4 cos (vt) + a5 cos (2vt) + b3 sin (vt) + b4 sin (2vt) , (24) η2 = a6 cos (vt) + a7 cos (2vt) + b5 sin (vt) + b6 sin (2vt) . By using the same method as above, we obtain ξ 2 = a3 + a5 cos (2vt) , (25) η2 = b6 sin (2vt) . Accordingly, the periodic solution of the Equation (14) with respect to e writes ξ = ξ 0 + cos (vt) e + [ a3 + a5 cos (2vt)] e2 , (26) η = b2 sin (vt) e + b6 sin (2vt) e2 , Universe 2020, 6, 110 15 of 25 where (ξ 0 , 0) is the coordinate of the collinear libration point. For Equation (15), suppose that the periodic solution in powers of parameter e takes the form ξ = ξ 1 e + ξ 2 e2 + ξ 3 e3 , (27) η = η1 e + η2 e 2 + η3 e 3 . Substituting Equation (27) into Equation (15) and comparing the terms with the same power of e. We obtain three sets of equations, of which the first two sets are identical to Equations (17) and (18) respectively, and the third set of equations has the form ξ̈ 3 − 2αω0 η̇3 = K1 ξ 3 + 2K2 ξ 1 ξ 2 + 2K3 η1 η2 + K4 ξ 13 + K5 ξ 1 η12 , (28) η̈3 + 2αω0 ξ̇ 3 = P1 η3 + P2 ξ 1 η2 + P2 ξ 2 η1 + P3 η13 + P4 ξ 12 η1 . Let the periodic solution of Equations (28) be ξ 3 = a8 cos (vt) + a9 cos (2vt) + a10 cos (3vt) + b7 sin (vt) + b8 sin (2vt) + b9 sin (3vt) , (29) η3 = a11 cos (vt) + a12 cos (2vt) + a13 cos (3vt) + b10 sin (vt) + b11 sin (2vt) + b12 sin (3vt) . Substituting Equations (23), (25) and (29) into Equation (28), and comparing the coefficients of sin ( jvt) and cos ( jvt) ( j = 1, 2, 3) respectively, then we have ξ 3 = a10 cos (3vt) , (30) η3 = b12 sin (3vt) , where the expressions of a10 and b12 can be found from Appendix B. Therefore, the periodic solution of Equation (15) in powers of parameter e becomes ξ = ξ 0 + cos (vt) e + [ a3 + a5 cos (2vt)] e2 + a10 cos (3vt) e3 , (31) η = b2 sin (vt) e + b6 sin (2vt) e2 + b12 sin (3vt) e3 , where the expression of coefficients, please refer to Appendix B again. Note that different values of κ not only determine the positions of the libration points, but also affect the shape of the second- and third-order periodic solutions through Equation (22). This is because the coefficients of Equation (22) are related to κ (see Appendix A for details). In order to better illustrate the approximate periodic solution obtained above, we show the phase portraits of the above approximate analytical periodic solutions near the libration point L3 under the four cases of κ = 0.5, 1, 2, and 4 respectively (see Figure 9a,d). When the value of κ is small, it is easy to find that the difference between the second- and the third-order approximate periodic solution is insignificant (see Figure 9a,b). Especially for the periodic solutions in Figure 9a, the second- and third-order approximate analytical periodic solutions on the two-dimensional plane around L3 agree very well, and they are almost indistinguishable visually. In this state, the second-order approximate periodic solution can be regarded as the alternative to the third-order case. If the values of other parameters remain unchanged, we can easily find from Figure 9c,d that the difference between the second- and third-order approximate periodic solutions become more and more significant as κ increases. Therefore, we can conclude according to Figure 9 that the smaller the value of κ, the better the degree of agreement between the second- and third-order approximate periodic solutions, and vice versa. Universe 2020, 6, 110 16 of 25 0.2 0.2 0 0 -0.2 -0.2 -1.3 -1.2 -1.1 -1.3 -1.2 -1.1 (a) (b) 0.4 0.5 0.2 0 0 -0.2 -0.5 -0.4 -1.25 -1.2 -1.15 -1.1 -1.3 -1.2 -1.1 -1 (c) (d) Figure 9. The approximate periodic solutions in two-dimensional plane for (a) κ = 0.5, (b) κ = 1, (c) κ = 2, and (d) κ = 4, respectively. 5.2. Expansion of Three-Dimensional Dynamic Equations Applying the translations u = u Li + ξ, v = η and w = ζ to Equation (7), then we have ∂Ω ξ̈ − 2αω0 η̇ = , ∂ξ ∂Ω η̈ + 2αω0 ξ̇ = , (32) ∂η ∂Ω ζ̈ = . ∂ζ Expanding the right hand side of Equation (32) to the second- and to the third-order respectively, then we obtain ξ̈ − 2αω0 η̇ = K1 ξ + K2 ξ 2 + K3 η 2 + K6 ζ 2 , η̈ + 2αω0 ξ̇ = P1 η + P2 ξη, (33) ζ̈ = Q1 ζ + Q2 ξζ, and ξ̈ − 2αω0 η̇ = K1 ξ + K2 ξ 2 + K3 η 2 + K6 ζ 2 + K7 ξ 3 + K8 ξη 2 + K9 ξζ 2 , η̈ + 2αω0 ξ̇ = P1 η + P2 ξη + P5 η 3 + P6 ξ 2 η + P7 ηζ 2 , (34) 3 2 2 ζ̈ = Q1 ζ + Q2 ξζ + Q3 ζ + Q4 ξ ζ + Q5 η ζ, Universe 2020, 6, 110 17 of 25 where the corresponding coefficients are displayed in the Appendix A. Three-Dimensional Periodic Solutions By means of successive approximations, we assume that the periodic solution of Equation (33) admits the following form ξ = ξ 1 e + ξ 2 e2 , η = η1 e + η2 e 2 , (35) 2 ζ = ζ1 e + ζ2 e . Substituting Equation (35) into Equation (33) and comparing these terms with the same power e, then we obtain ξ̈ 1 − 2αω0 η̇1 − K1 ξ 1 = 0, η̈1 + 2αω0 ξ̇ 1 − P1 η1 = 0, (36) ζ̈ 1 − Q1 ζ 1 = 0, and ξ̈ 2 − 2αω0 η̇2 − K1 ξ 2 − K2 ξ 12 − K3 η12 − K6 ζ 12 = 0, η̈2 + 2αω0 ξ̇ 2 − P1 η2 − P2 ξ 1 η1 = 0, (37) ζ̈ 2 − Q1 ζ 2 − Q2 ξ 1 ζ 1 = 0. For the third equation of Equation (36), the corresponding periodic solution takes the form ζ 1 = a cos (vt) + b sin (vt) , (38) and substitute Equation (38) into the third equation of Equation (36) yields   Q1 + v 2 ζ 1 = 0. (39) Note that if there is a non-zero solution to Equation (39), then Q1 = −v 2 < 0. So we choose a = 0 and b = 1, and one of the solutions is ζ 1 = sin (vt). Let the periodic solution of Equation (36) admits the following form ξ 1 = a14 cos (vt) + b13 sin (vt) , (40) η1 = a15 cos (vt) + b14 sin (vt) . Substituting Equation (40) into Equation (36) and comparing the coefficients of sin (vt) and cos (vt), then we obtain   K1 + v 2 a14 + 2αω0 vb14 = 0,   (41) 2αω0 va14 + P1 + v 2 b14 = 0, and   K1 + v 2 b13 − 2αω0 va15 = 0,   (42) − 2αω0 vb13 + P1 + v 2 a15 = 0. Universe 2020, 6, 110 18 of 25 Based on the results of Equations (20) and (21), then the periodic solution of Equation (36) writes ξ 1 = cos (vt) , η1 = b14 sin (vt) , (43) ζ 1 = sin (vt) . Similarly, let the periodic solution of Equation (37) has the following form ξ 2 = a16 + a17 cos (vt) + a18 cos (2vt) + b15 sin (vt) + b16 sin (2vt) , η2 = a19 cos (vt) + a20 cos (2vt) + b17 sin (vt) + b18 sin (2vt) , (44) ξ 2 = a21 cos (vt) + a22 cos (2vt) + b19 sin (vt) + b20 sin (2vt) . By using the same method as mentioned above, then we get ξ 2 = a16 + a18 cos (2vt) , η2 = b18 sin (2vt) , (45) ζ 2 = b20 sin (2vt) . Therefore, the periodic solution of Equation (33) with respect to e becomes ξ = ξ 0 + cos (vt) e + [ a16 + a18 cos (2vt)] e2 , η = b14 sin (vt) e + b18 sin (2vt) e2 , (46) 2 ζ = sin (vt) e + b20 sin (2vt) e , where the expressions of these coefficients can be obtained from Appendix B. As for the third-order approximate periodic solution of Equation (34), suppose that the periodic solution in the form of e can be written as ξ = ξ 1 e + ξ 2 e2 + ξ 3 e3 , η = η1 e + η2 e 2 + η3 e 3 , (47) 2 3 ζ = ζ1 e + ζ2 e + ζ3 e . Substituting Equation (47) into Equation (34), and comparing the terms containing the same power of e, then we obtain three set of equations, the equations corresponding to e and e2 are the same as Equations (36) and (37), respectively, and the equations corresponding to e3 are as follows ξ̈ 3 − 2αω0 η̇3 = K1 ξ 3 + 2K2 ξ 1 ξ 2 + 2K3 η1 η2 + 2K6 ζ 1 ζ 2 + K7 ξ 13 + K8 ξ 1 η12 + K9 ξ 1 ζ 12 , η̈3 + 2αω0 ξ̇ 3 = P1 η3 + P2 ξ 1 η2 + P2 ξ 2 η1 + P5 η13 + P6 ξ 12 η1 + P7 η1 ζ 12 , (48) ζ̈ 3 = Q1 ζ 3 + Q2 ξ 1 ζ 2 + Q2 ξ 2 ζ 1 + Q3 ζ 13 + Q4 ξ 12 ζ 1 + Q5 η12 ζ 1 . Writing the periodic solution of Equations (48) as ξ 3 = a23 cos (vt) + a24 cos (2vt) + a25 cos (3vt) + b21 sin (vt) + b22 sin (2vt) + b23 sin (3vt) , η3 = a26 cos (vt) + a27 cos (2vt) + a28 cos (3vt) + b24 sin (vt) + b25 sin (2vt) + b26 sin (3vt) , (49) ζ 3 = a29 cos (vt) + a30 cos (2vt) + a31 cos (3vt) + b27 sin (vt) + b28 sin (2vt) + b29 sin (3vt) . Substituting Equations (43), (45) and (49) into Equation (48), and comparing the coefficients of sin ( jvt) and cos ( jvt) ( j = 1, 2, 3), then we obtain Universe 2020, 6, 110 19 of 25 ξ 3 = a25 cos (3vt) , η3 = b26 sin (3vt) , (50) ζ 3 = b29 sin (3vt) . Therefore, the third-order approximate analytical periodic solution of Equation (7) in three-dimensional space can be written as ξ = ξ 0 + cos (vt) e + [ a16 + a18 cos (2vt)] e2 + a25 cos (3vt) e3 , η = b14 sin (vt) e + b18 sin (2vt) e2 + b26 sin (3vt) e3 , (51) 2 3 ζ = sin (vt) e + b20 sin (2vt) e + b29 sin (3vt) e . Note that one can find the specific expressions of the above coefficients from Appendix B. Based on the results of the above analysis, we show the three-dimensional periodic solution near the collinear libration point L3 for four different κ values of 0.5, 1, 2, and 4 in Figure 10. When κ = 0.5 and κ = 1, the second- and third-order three-dimensional periodic solutions agree very well, as shown in Figure 10a,b. When κ = 2, although Figure 10c shows that the second- and third-order periodic solutions are somewhat different, the degree of agreement is also good. However, Figure 10d (corresponding to κ = 4) indicates that the second- and third-order periodic solutions differ significantly. Therefore, it is easy to find that the smaller the value of κ, the better the approximation effect. (a) (b) (c) (d) Figure 10. The approximate periodic solutions in three-dimensional space for (a) κ = 0.5, (b) κ = 1, (c) κ = 2, and (d) κ = 4, respectively. Universe 2020, 6, 110 20 of 25 6. Conclusions For a restricted three-body problem consisting of a supergiant eclipsing binary system, we investigated the geometric structure of the corresponding zero-velocity surfaces and curves of the coplanar libration points and discussed the third body’s dynamic behavior under the influence of eight perturbation parameters through bifurcation diagrams. These perturbations include: Coriolis force α, centrifugal force β, oblateness coefficients α1 and α2 , random sums κ of the masses of the two primaries, mass ratio parameter µm , and the primaries’ radiation factors q1 and q2 . We found that any perturbation of Coriolis force and centrifugal force has a greater impact on the third body’s dynamic behavior, but their combined effect has a smaller effect on the third body. Similarly, the oblateness coefficient α2 has a more significant influence on the third body’s dynamic behavior than that of α1 . However, the degree of their combined effect on the third body’s dynamic behavior is minimal. Additionally, the parameter κ has a more extensive influence on the dynamic behavior of the third body, any tiny changes can cause apparent chaos phenomena. Furthermore, the smaller the mass ratio parameter, the more considerable the dynamic behavior of the third body, and vice versa. Besides, the influence of the radiation factor q1 on the three-body dynamics is more significant than that of q2 . The combined effect of q1 and q2 has a notable effect on the dynamics of the third body when the radiation factor q1 > 0.3434. In addition, we also applied the Lindstedt-Poincaré perturbation method to present the second- and third-order approximate analytical periodic solutions near the colinear libration points in the plane and space respectively, and illustrated the corresponding phase portraits with respect to κ = 0.5, 1, 2, and 4 around the libration point L3 . We found that the smaller the value of κ, the better the coincidence of the corresponding second-order and third-order periodic solutions. Besides, for the same κ value, the coincidence between the three-dimensional periodic solutions of different orders is better than that of the two-dimensional case. Author Contributions: Formal analysis, Y.W.; software, F.G.; writing-original draft, Y.W.; writing—review and editing, F.G. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the National Natural Science Foundation of China (NSFC) though grant No. 11672259, the China Scholarship Council through grant No. 201908320086. Acknowledgments: We are very grateful to the anonymous reviewers whose comments and suggestions helped improve and clarify this paper. Conflicts of Interest: The authors declare no conflict of interest. Appendix A ! ! 1 3α 1 3α K1 = ω02 ( β + κ − 1) + 2q1 κ (1 − µm ) 3 + 51 + 2q2 κµm 3 + 52 , d01 d01 d02 d02 ! ! 1 5α1 1 5α2 K2 = −3q1 κ (1 − µm ) + 6 V1 − 3q2 κµm + 6 V2 , d401 d01 d402 d02 ! ! 3 1 5α1 3 1 5α2 K3 = q1 κ (1 − µ m ) + 6 V1 + q2 κµm + 6 V2 , 2 d401 2d01 2 d402 2d02 ! ! 2 15α1 2 15α2 K4 = 2q1 κ (1 − µm ) + 7 + 2q2 κµm + 7 , d501 d01 d502 d02 ! ! 2 15α1 2 15α2 K5 = −3q1 κ (1 − µm ) + 7 − 3q2 κµm + 7 , d501 2d01 d502 2d02 ! ! 3 1 5α1 3 1 5α2 K6 = q1 κ (1 − µ m ) + 6 V1 + q2 κµm + 6 V2 , 2 d401 2d01 2 d402 2d02 Universe 2020, 6, 110 21 of 25 ! ! 2 15α 2 15α K7 = 2q1 κ (1 − µm ) 5 + 71 + 2q2 κµm 5 + 72 , d01 d01 d02 d02 ! ! 2 15α 2 15α K8 = −3q1 κ (1 − µm ) 5 + 71 − 3q2 κµm 5 + 72 , d01 2d01 d02 2d02 ! ! 2 15α 2 15α K9 = −3q1 κ (1 − µm ) 5 + 71 − 3q2 κµm 5 + 72 , d01 2d01 d02 2d02 ! ! 1 3α 1 3α P1 = ω02 ( β + κ − 1) − q1 κ (1 − µ m ) + 51 − q2 κµm + 52 , d301 2d01 d302 2d02 ! ! 1 5α 1 5α P2 = 3q1 κ (1 − µm ) + 61 V1 + 3q2 κµm + 62 V2 , d401 2d01 d402 2d02 ! ! 3 1 5α 3 1 5α2 P3 = q1 κ (1 − µm ) + 71 + q2 κµm + 7 , 2 d501 2d01 2 d502 2d02 ! ! 2 15α 2 15α2 P4 = −3q1 κ (1 − µm ) 5 + 7 1 − 3q2 κµm + 7 , d01 2d01 d502 2d02 ! ! 3 1 5α1 3 1 5α2 P5 = q1 κ (1 − µm ) 5 + 7 + q2 κµm + 7 , 2 d01 2d01 2 d502 2d02 ! ! 2 15α1 2 15α2 P6 = −3q1 κ (1 − µm ) + 7 − 3q2 κµm + 7 , d501 2d01 d502 2d02 ! ! 3 1 5α1 3 1 5α2 P7 = q1 κ (1 − µm ) 5 + 7 + q2 κµm + 7 , 2 d01 2d01 2 d502 2d02  V1 = sgn u Li + µm ,  V2 = sgn u Li + µm − 1 , d01 = u Li + µm , d02 = u Li + µm − 1 , ! ! 1 3α1 1 3α2 Q1 = ω02(κ − 1) − q1 κ (1 − µ m ) + 5 − q2 κµm + 5 , d301 2d01 d302 2d02 ! ! 1 5α 1 5α Q2 = 3q1 κ (1 − µm ) + 61 V1 + 3q2 κµm + 62 V2 , d401 2d01 d402 2d02 ! ! 3 1 5α1 3 1 5α2 Q3 = q1 κ (1 − µ m ) + 7 + q2 κµm + 7 , 2 d501 2d01 2 d502 2d02 ! ! 2 15α1 2 15α2 Q4 = −3q1 κ (1 − µm ) + 7 − 3q2 κµm + 7 , d501 2d01 d502 2d02 ! ! 3 1 5α1 3 1 5α2 Q5 = q1 κ (1 − µ m ) + 7 + q2 κµm + 7 . 2 d501 2d01 2 d502 2d02 Appendix B K2 + K3 b22 a3 = − , 2K1 Universe 2020, 6, 110 22 of 25 P1 + 4v 2 K3 b22 − K2 + 4αω0 vP2 b2   a5 = , 32v 4 + 8 K1 + P1 − 4α2 ω02 v 2 + 2K1 P1      1 1 1 1 P1 + 9v 2 −K2 a5 + K3 b2 b6 − K4 + K5 b22 + 3αω0 v P2 b6 + P2 a5 b2 − P3 b23 + P4 b2  4 4 2 2 a10 = , 81v 4 + 9 K1 + P1 − 4α2 ω02 v 2 + K1 P1  K2 + K3 b142 +K 6 a16 = − , 2K1 P1 + 4v 2 −K2 + K3 b14 2 + K + 4αω vP b   6 0 2 14 a18 = 4 2 2  2 , 32v + 8 K1 + P1 − 4α ω0 v + 2K1 P1   1 2  1 a25 = 4 2 2  2 P1 + 9v −K2 a18 + K3 b14 b18 + K6 b20 − K7 81v + 9 K1 + P1 − 4α ω0 v + K1 P1 4    1 2 + 1K 1 1 1 3 1 1 + K8 b14 9 − 6αω0 v − P2 b18 − P2 a18 b14 + P5 b14 − P6 b14 + P7 b14 , 4 4 2 2 4 4 4 2αω0 v b2 = − , P1 + v 2 − P2 b2 K1 + 4v 2 − 4αω0 v K3 b22 − K2   b6 = , 32v 4 + 8 K1 + P1 − 4α2 ω02 v 2 + 2K1 P1      1 1 1 1 1 1 K1 + 9v 2 − P2 b6 − P2 a5 b2 + P3 b23 − P4 b2 − 6αω0 v −K2 a5 + K3 b2 b6 − K4 + K5 b22  2 2 4 4 4 4 b12 = , 81v 4 + 9 K1 + P1 − 4α2 ω02 v 2 + K1 P1  2αω0 v b14 = − , P1 + v 2 − P2 b14 K1 + 4v 2 − 4αω0 v −K2 + K3 b14 2 +K   6 b18 = , 32v 4 + 8 K1 + P1 − 4α2 ω02 v 2 + 2K1 P1  Q2 b20 = − , 6v 2   1 2  1 1 1 3 b26 = 4 2 2  2 K1 + 9v − P2 b18 − P2 a18 b14 + P5 b14 81v + 9 K1 + P1 − 4α ω0 v + K1 P1 2 2 4    1 1 1 1 2 1 − P6 b14 + P7 b14 − 6αω0 v −K2 a18 + K3 b14 b18 + K6 b20 − K7 + K8 b14 + K9 , 4 4 4 4 4 1 1 1 1 1 2 − Q2 b20 − Q2 a18 + Q3 − Q4 + Q5 b14 b29 = 2 2 4 4 4 . 8v 2 References 1. Marchal, C. The Three-Body Problem; Elsevier Science: New York, NY, USA, 1990. 2. Ioka, K.; Chiba, T.; Tanaka, T.; Nakamura, T. Black hole binary formation in the expanding universe: Three body problem approximation. Phys. Rev. D 1998, 58, 063003. [CrossRef] 3. Nakamura, T.; Sasaki M.; Tanaka, T.; Thorne, K.S. Gravitational waves from coalescing black hole MACHO binaries. Astrophys. J. Lett. 1997, 487, L139. [CrossRef] 4. Zhou, T.Y.; Cao, W.G.; Xie, Y. Collinear solution to the three-body problem under a scalar-tensor gravity. Phys. Rev. D 2016, 93, 064065. [CrossRef] 5. Lindner, J.F.; Roseberry, M.I.; Shai, D.E.; Harmon, N.J.; Olaksen, K.D. Precession and chaos in the classical two-body problem in a spherical universe. Int. J. Bifurc. Chaos 2008, 18, 455–464. [CrossRef] 6. Gao, F.B.; Wang, R.F. Bifurcation analysis and periodic solutions of the HD 191408 system with triaxial and radiative perturbations. Universe 2020, 6, 35. [CrossRef] Universe 2020, 6, 110 23 of 25 7. Beletsky, V.V. Generalized restricted circular three-body problem as a model for dynamics of binary asteroids. Cosm. Res. 2007, 45, 408–416. [CrossRef] 8. Hou, X.Y.; Xin, X.S.; Feng, J.L. Forced motions around triangular libration points by solar radiation pressure in a binary asteroid system. Astrodynamics 2020, 4, 17–30. [CrossRef] 9. Braaten, E.; Kang, D.; Laha, R. Production of dark-matter bound states in the early universe by three-body recombination. J. High Energy Phys. 2018, 2018, 84. [CrossRef] 10. Trani, A.A.; Fujii, M.S.; Spera, M. The Keplerian three-body encounter. I. Insights on the origin of the S-stars and the G-objects in the galactic center. Astrophys. J. 2019, 875, 42. [CrossRef] 11. Muzzio, J.C.; Wachlin, F.C.; Carpintero, D.D. Regular and chaotic motion in a restricted three-body problem of astrophysical interest. Int. Astron. Union Colloq. 2000, 174, 281–285. [CrossRef] 12. Logoteta, D.; Bombaci, I. Constraints on microscopic and phenomenological equations of state of dense matter from GW170817. Universe 2019, 5, 204. [CrossRef] 13. Fahn, M.J.; Giesel, K.; Kobler, M. Dynamical properties of the Mukhanov-Sasaki Hamiltonian in the context of adiabatic vacua and the Lewis-Riesenfeld invariant. Universe 2019, 5, 170. [CrossRef] 14. Gao, F.B.; Zhang, W. A study on periodic solutions for the circular restricted three-body problem. Astron. J. 2014, 148, 116. [CrossRef] 15. Alshaery, A.A.; Abouelmagd E.I. Analysis of the spatial quantized three-body problem. Results Phys. 2020, 17, 103067. [CrossRef] 16. Zotos, E.E.; Chen, W.; Abouelmagd, E.I.; Han, H.T. Basins of convergence of equilibrium points in the restricted three-body problem with modified gravitational potential. Chaos Solitons Fractals 2020, 134, 109704. [CrossRef] 17. Alzahrani, F.; Abouelmagd, E.I.; Guirao, J.L.G.; Hobiny, A. On the libration collinear points in the restricted three-body problem. Open Phys. 2017, 15, 58–67. [CrossRef] 18. Abouelmagd, E.I.; Guirao, J.L.G.; Llibre, J. Periodic orbits for the perturbed planar circular restricted 3-body problem. Discret. Contin. Dyn. Syst.-Ser. B 2019, 24, 1007–1020. [CrossRef] 19. Selim, H.H.; Guirao, J.L.G.; Abouelmagd E.I. Libration points in the restricted three-body problem: Euler angles, existence and stability. Discret. Contin. Dyn. Syst.-Ser. B 2019, 12, 703–710. [CrossRef] 20. Abouelmagd E.I.; Mostafa A.; Guirao, J.L.G. A first order automated lie transform. Int. J. Bifurc. Chaos 2015, 25, 1540026.[CrossRef] 21. Abouelmagd, E.I.; Asiri, H.M.; Sharaf, M.A. The effect of oblateness in the perturbed restricted three-body problem. Meccanica 2013, 48, 2479–2490. [CrossRef] 22. Pathak, N.; Thomas, V.O.; Abouelmagd, E.I. The perturbed photogravitational restricted three-body problem: Analysis of resonant periodic orbits. Discret. Contin. Dyn. Syst.-Ser. B 2019, 12, 849–875. [CrossRef] 23. Niederman, L.; Pousse, A.; Robutel, P. On the co-orbital motion in the three-body problem: Existence of quasi-periodic horseshoe-shaped orbits. Commun. Math. Phys. 2020, 377, 551–612. [CrossRef] 24. Koon, W.S.; Lo, M.W.; Marsden J.E.; Ross, S.D. Dynamical Systems, the Three-Body Problem and Space Mission Design; Marsden Books: Wellington, New Zealand, 2011; ISBN 978-0-615-24095-4. 25. Qian, Y.J.; Yang, X.D.; Yang, L.Y.; Zhang, W. Approximate analytical methodology for the restricted three-body and four-body models based on polynomial series. Int. J. Aerosp. Eng. 2016, 2016, 9747289. [CrossRef] 26. Pathak, N.; Abouelmagd, E.I.; Thomas, V.O. On higher order resonant periodic orbits in the photo-gravitational planar restricted three-body problem with oblateness. J. Astronaut. Sci. 2019, 66, 475–505. [CrossRef] 27. Hénon, M. Generating Families in the Restricted Three-Body Problem; Springer: New York, NY, USA 1997. 28. Chenciner, A.; Montgomery, R. A remarkable periodic solution of the three-body problem in the case of equal masses. Ann. Math. 2000, 152, 881–901. [CrossRef] 29. Šuvakov, M.; Dmitrašinović, V. Three classes of Newtonian three-body planar periodic solutions. Phys. Rev. Lett. 2013, 110, 114301. [CrossRef] 30. Janković, M.R.; Dmitrašinović, V.; Šuvakov, M. A guide to hunting periodic three-body orbits with non-vanishing angular momentum. Comput. Phys. Commun. 2020, 250, 107052. [CrossRef] 31. Li, X.M.; Liao, S.J. More than six hundred new families of Newtonian periodic planar collisionless three-body orbits. Sci. China Phys. Mech. Astron. 2017, 60, 129511. [CrossRef] 32. Kalantonis, V.S. Numerical investigation for periodic solutions in the Hill three-body problem. Universe 2020, 6, 72. Universe 2020, 6, 110 24 of 25 33. Singh, J. Photogravitational restricted three-body problem with variable mass. Indian Indian J. Pure Appl. Math. 2003, 34, 335–341. 34. Singh, J.; Leke, O. Stability of the photogravitational restricted three-body problem with variable masses. Astrophys. Space Sci. 2010, 326, 305–314. [CrossRef] 35. Singh, J.; Leke, O. Equilibrium points and stability in the restricted three-body problem with oblateness and variable masses. Astrophys. Space Sci. 2012, 340, 27–41. [CrossRef] 36. Singh, J.; Leke, O. Effect of oblateness, perturbations, radiation and varying masses on the stability of equilibrium points in the restricted three-body problem. Astrophys. Space Sci. 2013, 344, 51–61. [CrossRef] 37. Singh, J.; Haruna, S. Equilibrium points and stability under effect of radiation and perturbing forces in the restricted problem of three oblate bodies. Astrophys. Space Sci. 2014, 349, 107–116.[CrossRef] 38. Singh, J.; Perdiou, A.E.; Gyegwe, J.M.; Kalantonis, V.S. Periodic orbits around the collinear equilibrium points for binary Sirius, Procyon, Luhman 16, α-Centuari and Luyten 726-8 systems: The spatial case. J. Phys. Commun. 2017, 1, 025008. [CrossRef] 39. Singh, J.; Perdiou, A.E.; Gyegwe, J.M.; Perdios, E.A. Periodic solutions around the collinear equilibrium points in the perturbed restricted three-body problem with triaxial and radiating primaries for binary HD 191408, Kruger 60 and HD 155876 systems. Appl. Math. Comput. 2018, 325, 358–374. [CrossRef] 40. Abdulraheem, A.R.; Singh, J. Combined effects of perturbations, radiation and oblateness on the periodic solutions in the restricted three-body problem. Astrophys. Space Sci. 2008, 317, 9–13.[CrossRef] 41. Abdulraheem, A.R.; Singh, J. Combined effects of perturbations, radiation, and oblateness on the stability of equilibrium points in the restricted three-body problem. Astron. J. 2006, 131, 1880. [CrossRef] 42. Singh, J.; Ishwar, B. Effect of perturbations on the location of equilibrium points in the restricted problem of three bodies with variable mass. Celestial Mech. 1984, 32, 297–305. [CrossRef] 43. Bhatnagar, K.B.; Chawla, J.M. A study of the Lagrangian points in the photogravitational restricted three-body problem. Indian J. Pure Appl. Math. 1979, 10, 1443–1451. 44. Szebehely, V. Stability of the points of equilibrium in the restricted problem. Astron. J. 1967, 72, 7–9. [CrossRef] 45. Subba Rao, P.V.; Sharma, R.K. Perturbations of the critical mass in the restricted three-body problem. Astrophys. Space Sci. 1976, 42, L17–L18. 46. Sharma, R.K. Perturbations of Lagrangian points in the restricted three-body problem. Indian J. Pure Appl. Math. 1975, 6, 1099–1102. 47. Sharma, R.K.; Subba Rao, P.V. Stationary solutions and their characteristic exponents in the restricted three-body problem when the more massive primary is an oblate spheroid. Celestial Mech. 1976, 13, 137–149. [CrossRef] 48. Abdullah, A.A.; Alhussain, Z.A.; Kellil, R. Locations and stability of the libration points in the CR3BP with perturbations. J. Math. Anal. 2017, 8, 131–144. 49. Abouelmagd, E.I.; Alzahrani, F.; Hobiny, A.; Guirao, J.L.G.; Alhothuali, M. Periodic solutions around the collinear libration points. J. Nonlinear Sci. Appl. 2016, 9, 1716–1727. [CrossRef] 50. Bekov, A.A. Particular solutions in the restricted collinear three body problem with variable masses. Soviet Astron. 1991, 35, 103–106. 51. Zotos, E.E. How does the oblateness coefficient influence the nature of solutions in the restricted three-body problem? Astrophys. Space Sci. 2015, 358, 33. [CrossRef] 52. Bhatnagar, K.B.; Hallan, P.P. Effect of perturbations in Coriolis and centrifugal forces on the stability of libration points in the restricted problem. Celestial Mech. 1978, 18, 105–112. [CrossRef] 53. Subbarao, P.V.; Sharma, R.K. A note on the stability of the triangular points of equilibrium in the restricted three-body problem. Astron. Astrophys. 1975, 43, 381–383. 54. Gyldén, H. Die bahnbewegungen in einem systeme von zwei Körpern in dem falle, dass die massen veränderungen unterworfen sind. Astron. Nachr. 1884, 109, 1–6. [CrossRef] 55. Bekov, A.A. Periodic solutions of the Gylden-Meshcherskii problem. Astron. Rep. 1993, 37, 651–654. 56. Bekov, A.A.; Momynov, S.B. Parametric solutions of the Gylden-Meshchersky problem. Int. J. Non-Linear Mech. 2019, 116, 195–199. [CrossRef] 57. Meshcherskii, I.V. Works on the Mechanics of Bodies of Variable Mass; GITTL: Moscow, Russia, 1952; p. 205. (In Russian) Universe 2020, 6, 110 25 of 25 58. Kalantonis, V.S.; Perdios, E.A.; Perdiou, A.E.; Vrahatis, M.N. Computing with certainty individual members of families of periodic orbits of a given period. Celestial Mech. Dyn. Astron. 2001, 80, 81–96. [CrossRef] 59. Hu, H.Y. Applied Nonlinear Dynamics; China Aeronautical Industry Press: Beijing, China, 2000. (In Chinese) c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).