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

Authors: Fabao Gao; Yongqing Wang
Journal: Advances in Astronomy
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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:(1)∑j=1∞Ajcos∑i=1kjiωit+Bjsin∑i=1kjiωit.The rate of change of the solution with respect to timet in equation (1) appears as a combination of many incommensurable frequencies, ωii=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. The 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. The 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. They 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. Thirteen 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. They 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 second- and 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. The 3D multiple time scales solution has the following form:(2)u=∑l=12εlult,εt,where u=x,y,z, ul=x1,y1,z1, l=1,2, and ε is a small, dimensionless parameter.


## 2. 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:(3)x¨−2y˙=x−Ux,y¨+2x˙=y−Uy,z¨=−UZ,where the potential function(4)Ux,y,z=−1−μx+μ2+y2+z21/2−μx−1+μ2+y2+z21/2.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. Therefore, we expect that the periodic solution which is beneficial to the actual mission can be realized by “construction.”Figure 1
(a–c) Bifurcation diagrams of mass ratio parameters inµ-x-y, µ-x-z, and µ-y-z frames, respectively.
(a)(b)(c)Next, we first deal with the complex potential functionUx,y,z by introducing Legendre polynomials Pn, which satisfies Rodrigues formula:(5)Pnx=12nn!dndxnx2−1n,n=0,1,2,3,….Then, we have [22](6)1x−A2+y−B2+z−C21/2=1D∑n=0∞ρDnPnAx+By+CzDρ,where D2=A2+B2+C2 and ρ=x2+y2+z2.Substituting equations (4) and (6) into equation (3), it yields(7)x¨−2y˙−1+2c2x=∂∂x∑n≥3cnρnPnxρ,y¨+2x˙+c2−1y=∂∂y∑n≥3cnρnPnxρ,z¨+c2z=∂∂z∑n≥3cnρnPnxρ,where cn=1/γL3±1nμ+−1n1−μγLn+1/1∓γLn+1, for libration points L1 and L2 and cn=1/γL31−μ+μγLn+1/1+γLn+1, for libration point L3, γL=n12/3, where n1 represents the angular velocity of relative motion between the two primaries.Then, equation (7), up to the third approximation, can be written in the following form:(8)x¨−2y˙−1+2c2x=32c32x2−y2−z2+2c42x2−3y2−3z2,y¨+2x˙+c2−1y=−3c3xy−32c4y4x2−y2−z2,z¨+c2z=−3c3xz−32c4z4x2−y2−z2.


## 3. Construction of Approximate 3D Multiple Time Scales Solution

Consider the small amplitude of nonhomogeneous terms of equation (8). Then,(9)x¨−2y˙−1+2c2x=32εc32x2−y2−z2+2εc42x2−3y2−3z2,y¨+2x˙+c2−1y=−3εc3xy−32εc4y4x2−y2−z2,z¨+c2z=−3εc3xz−32εc4z4x2−y2−z2.We write the solution of equation (9) as follows:(10)x=∑l=12εlxlT0,T1,y=∑l=12εlylT0,T1,z=∑l=12εlzlT0,T1,where T0=tandT1=εt.Then, time derivatives become(11)ddt=∂∂T0dT0dt+∂∂T1dT1dt+⋯=D0+εD1+⋯,(12)d2dt2=D0+εD1+⋯2=D02+2εD0D1+⋯,where Dk=∂/∂Tk,k=0,1,….Substituting equations (10)–(12) into equation (9), we find(13a)D02x1−2D0y1−1+2c2x1ε+D02x2+2D0D1x1−2D0y2−2D1y1−1+2c2x2ε2+2D0D1x2−D1y2ε3=32c32x12−y12−z12ε3+2c42x12−3y12−3z12x1+3c32x1x2−y1y2−z1z2ε4+Οε5,(13b)D02y1+2D0x1+c2−1y1ε+D02y2+2D0D1y1+2D0x2+2D1x1+c2−1y2ε2+2D0D1y2+D1x2ε3=−3c3x1y1ε3−3c3x1y2+x2y1+12c44x12−y12−z12y1ε4+Οε5,(13c)D02z1+c2z1ε+D02z2+2D0D1z1+c2z2ε2+2D0D1z2ε3=−3c3x1z1ε3−3c3x1z2+x2z1+12c44x12−y12−z12z1ε4+Οε5.Now, from equations (13a)–(13c), equating the coefficients of ε, ε2, and ε3 to zero, it leads to the following three sets of equations.Orderε:(14)D02x1−2D0y1−1+2c2x1=0,D02y1+2D0x1+c2−1y1=0,D02z1+c2z1=0.Orderε2:(15)D02x2+2D0D1x1−2D0y2−2D1y1−1+2c2x2=0,D02y2+2D0D1y1+2D0x2+2D1x1+c2−1y2=0,D02z2+2D0D1z1+c2z2=0.Orderε3:(16)2D0D1x2−2D1y2=32c32x12−y12−z12,2D0D1y2+2D1x2=−3c3x1y1,2D0D1z2=−3c3x1z1.Consider that the last equation of (14) is decoupled from the first two equations of (14). We employ the following transformations:(17)x1=ω1,ω1˙=ω2,y1=ω3,ω3˙=ω4.Then, the first two equations of (14) can be represented by the following equivalent form:(18)ω1˙=ω2,ω2˙=1+2c2ω1+2ω4,ω3˙=ω4,ω4˙=−2ω2+1−c2ω3.Thus, the coefficient matrix of equation (18) can be denoted as(19)01001+2c200200010−21−c20.Therefore, the characteristic equation can be obtained as(20)λ^4+2−c2λ^2+1+2c21−c2=0,where(21)λ^1=12−4+2c2+29c22−8c21/21/2,λ^2=−12−4+2c2+29c22−8c21/21/2,λ^3=12−4+2c2−29c22−8c21/21/2,λ^4=−12−4+2c2−29c22−8c21/21/2.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:(22)x1=I1eλ^1t+I2e−λ^2t+I3cosλt+I4sinλt,y1=−k1I1eλ^1t+k1I2e−λ^2t−k2I3sinλt+k2I4cosλt,z1=J1cosωt+J2sinωt,where coefficients I1,…,I4,k1,k2,J1 and J2 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 thatI1=I2=0, which implies that(23)x1=I3cosλt+I4sinλt,y1=−k2I3sinλt+k2I4cosλt,z1=J1cosωt+J2sinωt,i.e.,(24)x1=−Axcosλt+ϕ,y1=k2Axsinλt+ϕ,z1=Azsinωt+ψ,where(25)Ax=I32+I42,AZ=J12+J22,k2=1+2c2+λ22λ=2λλ2−c2+1,cosϕ=−I3Ax,sinϕ=I4Ax,cosψ=−J2AZ,sinψ=J1AZ.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). The 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. Therefore, 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.Figure 2
Time history diagrams with respect to (a)x and (b) y, respectively, when μ=0.0123.
(a)(b)Figure 3
Quasiperiodic orbit of Lissajous type inx-y-t space for μ=0.0123.Figure 4
Quasiperiodic orbit of Lissajous type in the (a)x-y plane and (b) x-y-z space, respectively, when μ=0.0123.
(a)(b)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. Therefore, without loss of generality, we consider here the case ofλ=ω. Then, equation (24) can be represented as(26)x1=A1T1eiλT0+A¯1T1e−iλT0,y1=ik2A1T1eiλT0−ik2A¯1T1e−iλT0,z1=A2T1eiλT0+A¯2T1e−iλT0.Substituting equation (26) into equation (15), it yields(27)D02x2−2D0y2−1+2c2x2+2iλD1A1eiλT0−2iλD1A¯1e−iλT0−2ik2D1A1eiλT0+2ik2D1A¯1e−iλT0=0,D02y2+2D0x2+c2−1y2−2λk2D1A1eiλT0−2λk2D1A¯1e−iλT0+2D1A1eiλT0+2D1A¯1e−iλT0=0,D02z2+c2z2+2iλD1A2eiλT0−2iλD1A¯2e−iλT0=0.Let the secular term be equal to zero, and we find that(28)D1A1λ−k2=0,D1A1k2−1=0,D1A2=0.Hence,(29)A1=ρ12,A2=ρ22,where ρ1 and ρ2 are constants.Substituting equation (29) into equation (26), and set constant k2=δ, we obtain(30)x1=ρ1cosλt,y1=−δρ1sinλt,z1=ρ2cosλt.Then, equation (27) can be expressed as(31)D02x2−2D0y2−1+2c2x2=0,D02y2+2D0x2+c2−1y2=0,D02z2+c2z2=0.Similarly, solutions to equation (31) can be denoted as(32)x2=B1T1eiλT0+B¯1T1e−iλT0,y2=iδB1T1eiλT0−iδB¯1T1e−iλT0,z2=B2T1eiλT0+B¯2T1e−iλT0.Substituting equations (26) and (32) into equation (16), it yields(33)2iλ−δD1B1eiλT0+2iδ−λD1B¯1e−iλT0=34c32+δ2ρ12−ρ22cos2λT0+2−δ2ρ12−ρ22,21−λδD1B1eiλT0+21−λδD1B¯1e−iλT0=−3c3iδA12ei2λT0+3c3iδA¯12e−i2λT0,2iλD1B2eiλT0−2iλD1B¯2e−iλT0=−3c3A1A2ei2λT0−3c3A¯1A¯2e−i2λT0−3c3A1A¯2+A¯1A2.We define(34)B1T1=12α1T1eiβ1T1,B2T1=12α2T1eiβ2T1.Thus,(35)D1B1=12α1′eiβ1+12α1eiβ1iβ1′,D1B¯1=12α1′e−iβ1+12α1e−iβ1−iβ1′,D1B2=12α2′eiβ2+12α1eiβ2iβ2′,D1B¯2=12α2′e−iβ2+12α2eiβ2−iβ2′.Substituting equations (35) and (30) into equation (33), we get(36)2λ−δα1′sinλT0+β1+α1β1′cosλT0+β1=34c32+δ2ρ12−ρ22cos2λT0+2−δ2ρ12−ρ22,21−λδα1′cosλT0+β1−α1β1′sinλT0+β1=32c3δρ12sin2λT0,−2λα2′sinλT0+β2+α2β2′cosλT0+β2=−32c3δρ1ρ2cos2λT0+1.Then, it can be followed from equation (36) that(37)α1′T1=38δ−λc32+δ2ρ12−ρ22cos2λt+2−δ2ρ12−ρ22sinλt+β1T1+341−λδc3δρ12sin2λtcosλt+β1T1,α1β1′T1=38δ−λc32+δ2ρ12−ρ22cos2λt+2−δ2ρ12−ρ22cosλt+β1T1−341−λδc3δρ12sin2λtsinλt+β1T1.Substituting equation (34) into equation (32), we arrive at(38)x1=α1T1cosλt+β1T1,y1=−δα1T1sinλt+β1T1,z1=α2T1cosλt+β2T1.Combining equations (10), (30), and (38), an approximate 3D multiple time scales solution to equation (9) can be represented as(39)x=ερ1cosλt+ε2cosλt+β1εt,y=−εδρ1sinλt−ε2δα1εtsinλt+β1εt,z=ερ2cosλt+ε2α2εtcosλt+β2εt.where α1 and β1 satisfy equation (37) and α2 and β2 satisfy the third equation of (36).


## 4. 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. The 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. The 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. Therefore, 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.



---

*Source: 8868137-2020-11-11.xml*