Approximate Analytical Periodic Solutions to the Restricted Three-Body Problem with Perturbation, Oblateness, Radiation and Varying Mass

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.


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][10][11][12][13] for more information).
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.

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 where r 2 i = (x − x i ) 2 + y 2 + z 2 (i = 1, 2) denotes the distance between the ith primary and the third body, and the two primaries are placed on (x 1 , 0, 0) and (x 2 , 0, 0), respectively. For the distance r between the primaries, represents the mass of the ith primary, f is the gravitational constant. In addition, the parameters q i and A i (t) stand for the radiation factor and the varying oblateness of the ith primary, respectively. The Coriolis force α and centrifugal force β satisfy α = 1 + ε, β = 1 + ε , here |ε| and |ε | 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 where κ is an arbitrary dimensionless constant of a particular integral rµ(t) = κC 2 in Gyldén-Meshcherskii problem (see References [54][55][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, A 2 (t) = 0. Then Equation (2) is simplified to ω 2 = 1 + 3A 1 (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]) where d i is the distance between the third body and the ith primary, R(t) = φt 2 + 2ϕt + γ, φ, ϕ and γ are constants, the unified Meshcherskii law where µ 0 , µ 10 , µ 20 are constants, and the particular solutions of Gyldén-Meshcherskii problem, where d 12 is the distance between the two primaries at initial time. Applying transformation A i (t) = 2) refers to the oblateness of the ith primary at initial time, d E i and d P i 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 = κ, d 12 = 1. Therefore, we obtain It follows from Equation (5) that ω 2 0 = 1 + 3(α 1 + α 2 )/2, α i 1 (i = 1, 2). Then Equations (1) are reduced to the following formü with the potential function where are the distances from the third body to the two primaries, respectively.

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 where C is the Jacobi constant andu 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 q 1 = q 2 = 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 L 1 , L 2 , · · · , L 5 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 L 1 ), and the third body can pass from one primary to the other through L 1 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" L 2 , 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 L 3 . 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, the symmetrical libration points L 4 and L 5 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 L 1 (0.018995, 0), L 2 (1.1989, 0), L 3 (−1.19261, 0), L 4 (0.01424987, 0.851953), L 5 (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].

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 u 0 = (−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).  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 q 1 and q 2 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 q 1 or q 2 has almost the same influence on the third body's movement in u and v directions. However, the influence of q 1 on the third body's amplitude is greater than that of q 2 , especially when q 1 < 0.57 or 0.67 < q 1 < 0.85. In contrast, the impact of q 2 is much more moderate. Figure 8c shows that when q 1 ≤ 0.0202, regardless of how q 2 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 q 1 ∈ [0.0202, 0.3434], the third body's amplitude has a significant change, even if q 2 is very small. When the value of q 1 is greater than 0.3434, the third body exhibits approximately symmetrical motion in a larger area under the joint action of q 1 and q 2 .    (c) Figure 8. (a-c) Bifurcation diagrams of radiation factor of the ith (i = 1, 2) primary in the frames q 1 uv, q 2 uv and q 1 q 2 u, respectively, when κ = 1, µ m = 0.48785, α 1 = 0.024, α 2 = 0.02, α = 1.001, β = 1.002.

Expansion of Two-Dimensional Dynamic Equations
The third body's equations of motion in the u-v plane arë where the potential function where d a = (u + µ m ) 2 + v 2 and d b = (u + µ m − 1) 2 + v 2 are the distances from the third body to two primaries respectively. Then the right hand side of Equation (10) can be written as Now we study the periodic solutions of the third body around the collinear libration points L i (i = 1, 2, 3), and u L i is the specific coordinate corresponding to L i . By letting u = u L i + ξ and v = η, then Equation (10) are reduced to the following form Expanding the right hand side of Equation (13) to the second-order yields ξ − 2αω 0η = K 1 ξ + K 2 ξ 2 + K 3 η 2 , η + 2αω 0ξ = P 1 η + P 2 ξη, (14) and to the third-order in the same way, then we obtain where corresponding coefficients can be found in Appendix A.
Substituting Equation (16) into Equation (14) and comparing the terms with the same power of e, we have 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 with period T = 2π/ .
Substituting Equation (19) into Equation (17) and comparing the coefficients of sin ( t) and cos ( t), then we get and According to the homogeneous linear Equation (20) and (21), in order to get a nonzero solution, the corresponding determinant satisfies then we set a 1 = 1, b 1 = 0, and substitute them into Equations (20) and (21) respectively, such that a 2 = 0 and b 2 = −2αω 0 P 1 + 2 , so the periodic solution of Equation (17) is Similarly, we write the periodic solution of Equation (18) as By using the same method as above, we obtain ξ 2 = a 3 + a 5 cos (2 t) , Accordingly, the periodic solution of the Equation (14) with respect to e writes ξ = ξ 0 + cos ( t) e + [a 3 + a 5 cos (2 t)] e 2 , 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 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 Let the periodic solution of Equations (28) be Substituting Equations (23), (25) and (29) into Equation (28), and comparing the coefficients of sin (j t) and cos (j t) (j = 1, 2, 3) respectively, then we have ξ 3 = a 10 cos (3 t) , where the expressions of a 10 and b 12 can be found from Appendix B.
Therefore, the periodic solution of Equation (15) in powers of parameter e becomes ξ = ξ 0 + cos ( t) e + [a 3 + a 5 cos (2 t)] e 2 + a 10 cos (3 t) e 3 , 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 L 3 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 L 3 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.

Three-Dimensional Periodic Solutions
By means of successive approximations, we assume that the periodic solution of Equation (33) admits the following form Substituting Equation (35) into Equation (33) and comparing these terms with the same power e, then we obtainξ For the third equation of Equation (36), the corresponding periodic solution takes the form and substitute Equation (38) into the third equation of Equation (36) yields Note that if there is a non-zero solution to Equation (39), then Q 1 = − 2 < 0. So we choose a = 0 and b = 1, and one of the solutions is ζ 1 = sin ( t).
Let the periodic solution of Equation (36) admits the following form Substituting Equation (40) into Equation (36) and comparing the coefficients of sin ( t) and cos ( t), then we obtain and Based on the results of Equations (20) and (21), then the periodic solution of Equation (36) writes Similarly, let the periodic solution of Equation (37) has the following form By using the same method as mentioned above, then we get ξ 2 = a 16 + a 18 cos (2 t) , Therefore, the periodic solution of Equation (33) with respect to e becomes ξ = ξ 0 + cos ( t) e + [a 16 + a 18 cos (2 t)] e 2 , where the expressions of these coefficients can be obtained from Appendix B.
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 e 2 are the same as Equations (36) and (37), respectively, and the equations corresponding to e 3 are as follows Writing the periodic solution of Equations (48) as Substituting Equations (43), (45) and (49) into Equation (48), and comparing the coefficients of sin (j t) and cos (j t) (j = 1, 2, 3), then we obtain ξ 3 = a 25 cos (3 t) , Therefore, the third-order approximate analytical periodic solution of Equation (7) in three-dimensional space can be written as ξ = ξ 0 + cos ( t) e + [a 16 + a 18 cos (2 t)] e 2 + a 25 cos (3 t) e 3 , 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 L 3 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.

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 q 1 and q 2 . 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 q 1 on the three-body dynamics is more significant than that of q 2 . The combined effect of q 1 and q 2 has a notable effect on the dynamics of the third body when the radiation factor q 1 > 0.3434.
In addition, we also applied the Lindstedt-Poincaré perturbation method to present the secondand 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 L 3 . 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. , , V 1 = sgn u L i + µ m , .