Approximate Analytical Three-Dimensional Multiple Time Scales Solution to a Circular Restricted Three-Body Problem

We illustrate the chaotic nature of the circular restricted three-body problem from the perspective of the bifurcation diagram with respect to the mass ratio parameter. Moreover, it is shown that when the frequency ratio in different directions of the classical problem is irrational, it has the quasiperiodic characteristics. In addition, a three-dimensional approximate solution to this problem under two time scales is proposed by using the multiple time scales method.


Introduction
As early as the 19th century, plenty of mathematicians, such as Dirichlet and Weierstrass, et al. [1], expected to obtain a series solution of the three-body problem in the following form: e rate of change of the solution with respect to time t in equation (1) appears as a combination of many incommensurable frequencies, ω i i � 1, 2, . . . , k, often called quasiperiodic solution.
Almost half a century ago, Farquhar and Kamel [2] proposed an approximation to construct the periodic orbits described above. A few years later, Richardson [3,4] developed third-order analytic solutions to collinear libration points (CLP) for a class of circular restricted three-body problem (CRTBP) based on a method similar to Lindstedt-Poincaré and successive approximation. e constructed solutions are the basis of determining halo orbits around these points, and they can provide approximate initial values of halo orbits. Furthermore, the solutions can also be developed to investigate spacecraft formation flying. e more accurate the approximate analytical solutions, the less fuel in the system will be consumed. Ruijgrok [5] solved a planar three-body problem for equal mass particles under a particular class of three-body forces. Zagouras and Markellos [6] studied a periodic spatial solution to Hill's problem, which is the limiting case of a restricted three-body problem (R3BP). Papadakis et al. [7] considered the periodic orbits of the three-body problem generated by CLP and provided the second-order approximate analytic expressions of such orbits. Gómez et al. [8] computed quasihalo orbits in CRTBP semianalytically through the Lindstedt-Poincaré technique. Ten years ago, Lu and Zhao [9] proposed an improved and more precise third-order analytic solution than Richardson's classical analytic solution. Gidea and Deppe [10] numerically studied the influence of small perturbations on the dynamics of an infinitesimal third body. ey investigated the chaotic orbits in an RTBP, which also indicated the complexity and difficulty of obtaining approximate analytic solutions to the problem.
In recent years,Šuvakov and Dmitrašinovic [11] presented some creative results for the periodic orbits of the three-body problem.
irteen new, zero-angular-momentum, distinct equal mass, planar collisionless periodic orbits are displayed in three new categories. When the two primaries are both oblate spheres, Mittal et al. [12] used the predictor-corrector algorithm to construct periodic orbits. ey showed the corresponding periodic orbits under different energy constants, mass ratios, and oblateness factors of the two primaries. Based on photogravitational planar RTBP with oblateness, Pathak et al. [13] studied the seventh-, ninth-, and eleventh-order internal resonance periodic orbits of the Sun-Earth system by using Runge-Kutta-Gill method. For a generalized photogravitational RTBP model, the two primaries are oblate spheroid and are under the gravity of an asteroid belt. Abouelmagd et al. [14] found a secular solution that can be reduced to a periodic one near the triangular libration points when the mass ratio is equal to the critical mass value. Similar results can also be found in [15]. Selim et al. [16] analyzed the stable motion solutions of long-period orbits and short-period orbits when the primary is a triaxial rigid body and the Euler angle of rotation satisfies certain conditions. Besides, for two types of the perturbed RTBP, Gao et al. [17,18] utilized the Lindstedt-Poincaré perturbation method to give an approximate analytical solution to the periodic orbits near the CLP.
Considering that the dynamical equation of the third body in CR3BP is a time-varying high-dimensional nonlinear system in the inertial frame, it can hardly be solved by using analytical approaches. However, this case will be better in the rotating coordinate system, for the aforementioned governing equations are time-independent, and there is an integrable motion. Using the method of multiple time scales, Nayfeh [19,20] studied the 3 : 1 and 2 : 1 small-amplitude resonances near the triangular libration points in the plane when the potential function was expanded to the secondand third-order terms of a small parameter, respectively.
Based on the number of terms expanded by the potential function of the problem, the existing time-scale solutions mainly include the three-dimensional single time-scale solution when the potential function expands to the third-order, the three-dimensional (3D) double time-scale solution when the potential function expands to the second-order, and the planar approximate analytic solution. Accordingly, when the third body is also considered to vibrate with a small amplitude in space, the multiple time scales method [21] will be used to construct the approximate 3D CR3BP solution in this paper. e 3D multiple time scales solution has the following form: where u � (x, y, z), u l � (x 1 , y 1 , z 1 ), l � 1, 2, and ε is a small, dimensionless parameter.

Dynamic Equations of Classical CRTBP
In this section, the dynamic equations of classical CRTBP in three-dimensional space will be described.
In a rotating framework, it is well known that the governing equations of the classical CRTBP can be characterized by the following set of differential equations: where the potential function For more details on the establishment of equation (3), please refer to [22].
It is easy to find that μ in equation (4) is the only parameter in the classical CRTBP. In this paper, we define it as the ratio of the mass of one primary to the mass sum of the two primaries, so it is called the mass ratio parameter, and its value range is between 0 and 1. Under the initial conditions [−0.001, 0.2, 0.1, −0.3, 0.87, 0.01], we divide the range of μ into 1000 equal parts, and then, respectively, return its corresponding coordinates of three different planes (x-y, x-z, y-z) for each value of μ, and the bifurcation diagram with respect to the parameter μ based on ode45 algorithm shows that any small change of μ will lead to the third body's apparent chaotic dynamic behavior (see Figure 1). We find that any value of μ corresponds to an infinite number of points rather than a single point, which means that the probability of obtaining periodic solutions of the three-body problem by numerical simulation is almost zero. Moreover, even if there are some results that are seem to be periodic, part of the existing conclusions based on numerical simulation will no longer hold when the step size continues to shrink. When the iteration accuracy is gradually improved, it is difficult for the orbit of the three-body system with chaotic characteristics to come back after iteration from a particular initial point, and researchers do not even know what will happen when the accuracy is less than 10 −16 . erefore, we expect that the periodic solution which is beneficial to the actual mission can be realized by "construction." Next, we first deal with the complex potential function U(x, y, z) by introducing Legendre polynomials P n , which satisfies Rodrigues formula: en, we have [22] 1 where Substituting equations (4) and (6) into equation (3), it yields where , for libration points L 1 and L 2 and where n 1 represents the angular velocity of relative motion between the two primaries.
en, equation (7), up to the third approximation, can be written in the following form:

Construction of Approximate 3D Multiple Time Scales Solution
Consider the small amplitude of nonhomogeneous terms of equation (8). en, We write the solution of equation (9) as follows: where T 0 � t and T 1 � εt.

Advances in Astronomy
Now, from equations (13a)-(13c), equating the coefficients of ε, ε 2 , and ε 3 to zero, it leads to the following three sets of equations.
Order ε: Order ε 2 : Order ε 3 : Consider that the last equation of (14) is decoupled from the first two equations of (14). We employ the following transformations: en, the first two equations of (14) can be represented by the following equivalent form: us, the coefficient matrix of equation (18) can be denoted as erefore, the characteristic equation can be obtained as where It is noticed that λ 1 and λ 2 are two equal and opposite real roots, and λ 3 and λ 4 are a pair of pure imaginary roots. Hence, the general solution of equation (14) assumes the following form: x 1 � I 1 e λ 1 t + I 2 e − λ 2 t + I 3 cos(λt) + I 4 sin(λt), y 1 � −k 1 I 1 e λ 1 t + k 1 I 2 e − λ 2 t − k 2 I 3 sin(λt) + k 2 I 4 cos(λt), where coefficients I 1 , . . . , I 4 , k 1 , k 2 , J 1 and J 2 are all constants and λ is the mould of λ 3 and λ 4 .
Since the motion of the third body may be unbounded due to the influence of exponential terms, appropriate initial conditions are selected such that I 1 � I 2 � 0, which implies that x 1 � I 3 cos(λt) + I 4 sin(λt), y 1 � −k 2 I 3 sin(λt) + k 2 I 4 cos(λt), i.e., Advances in Astronomy where If the value of (λ/ω) in equation (24) is an irrational number, according to the time history diagrams in Figure 2, this linearization motion of the system appears to make the periodic motion in the x and y components, respectively. However, it is not periodic but Lissajous-type quasiperiodic in the x-y plane, and this fact can be verified by the integral curves with respect to x and y (see Figure 3) and its projections in the x-y plane, as well as the portraits in the frame of x-y-z (see Figure 4). e significance of Figures 2-4 is that if the frequency of the third body in all directions is not commensurable, it will exhibit a quasiperiodic motion. erefore, the value of (λ/ω) is a rational number that performs a critical role in the motion law of the system. In general, if the value of (λ/ω) is an irrational number, the time history diagram can only serve as a reference but a standard to judge if the system makes the periodic motion or not. On the contrary, the motion law of the system can be judged by phase portraits and time history diagrams.
In addition, note that the halo-type periodic orbits play particular importance in the actual mission ACE, ISEE-3/ ICE, MAP, SOHO, Genesis, etc., that is, the in-plane and out-of-plane motion makes those characteristic frequencies equal within the sufficiently large motion region. erefore, without loss of generality, we consider here the case of λ � ω.

Conclusions
In this paper, we first show the chaotic characteristics of the classical CRTBP through the bifurcation diagram concerning the mass parameter μ, which indicates that it is challenging to find the periodic orbit by adjusting the initial conditions of the problem. e astronomer Poincaré has also shown that the probability of the event is zero. Secondly, we also show that even though sometimes the third body may exhibit corresponding periodic motions in different directions in space, the incommensurability of their frequencies will cause the third body to have quasiperiodic characteristics.
When the third body moves with a small amplitude, we construct a kind of approximate 3D multiscale solution by decomposing the conventional time contained in the solution of the CRTBP into the coupling effect of multiple time scales. e shortcoming of this type of solution is that, like other known approximate solutions (such as those based on the Lindstedt-Poincaré method), the result of error analysis is not satisfactory because of the chaotic nature of the problem. In theory, we need more infinite terms in solution series and time scales to make it approach the exact solution of the problem, but this is worthless, and the analytical calculation under the current technological level is almost impossible to achieve. erefore, we only need the series approximation of two or three terms, so that the multiple time scales solution can be used as an initial guess or approximation of the low-energy transfer orbital mission (such as Genesis discovery mission), and then generate spacecraft orbit through some differential correction procedures. In addition, the classic ISEE-3 mission actually used a similarly constructed approximate solution as the initial solution for orbit design.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.