Approximate Analytical Methodology for the Restricted Three-Body and Four-Body Models Based on Polynomial Series

The restricted three-body problem (R3BP) and restricted four-body problem (R4BP) are modeled based on the rotating frame.The conservative autonomous system for the R3BP and nonautonomous system with period parametric resonance due to the fourth body are derived. From the vibrational point of view, the methodology of polynomial series is proposed to solve for these problems analytically. By introducing the polynomial series relations among the three directions of motion, the three-degree-of-freedom coupled equations are transferred into one degree-of-freedom containing the full dynamics of the original autonomous system for the R3BP. As for the R4BP case, the methodology of polynomial series combined with the iterative approach is proposed. During the iterative approach, the nonautonomous system can be treated as pseudoautonomous equation and the final polynomial series relations and one-degree-of-freedom system can be derived iteratively.


Introduction
In our solar system, the planets move around the Sun, while the moons move around their host planets, which in turn also move around the Sun.They all follow the two-body theory or the perturbed two-body theory which has well-known closed form solution given in terms of elementary functions.However, to study the motion of artificial satellites in the solar system, the problem is too complicated to obtain similar types of solutions.Compared to the mass of either the Sun, the planet, or the moon, the mass of the artificial satellite is always infinite, which means that the gravitational influence of the artificial satellite on the primaries can be omitted.According to the different mission purposes, the system can be modeled correspondingly as the restricted three-body problem (R3BP) or the restricted four-body problem (R4BP).For example, if the motion of the satellite is limited in the planet-moonmoon system, such as the Saturn-Tethys-Telesto-spacecraft, Saturn-Tethys-Calypso-spacecraft, and Saturn-Dione-Helenspacecraft systems, they are typical R4BP systems.If the motion of the satellite is limited in the planet-moon system while the Sun's gravitational influence cannot be neglected, such as ARTEMIS mission [1], the system needs to be modeled as the R4BP as well.Meanwhile, for the mission exploring the Sun-planet system, such as GALA and PLANCK [2], the R3BP model is suitable when we consider the Sun and the barycenter of the planet and its moon as two dominant masses move around their center of mass along either circular or elliptic orbits.
At first glance, the difficulties of the R3BP and R4BP are not obvious, especially when considering that the two-body problem can be analytically solved.In the past, many physicists, astronomers, and mathematicians attempted unsuccessfully to find closed form solutions to the three-body problem.Such solutions do not exist because motions of the bodies are in general unpredictable, which makes the problem one of the most challenging problems in the history of science.
The analytical techniques developed for studying the dynamics around the CR3BP and the R4BP models are similar.Although the exact solutions for R3BP and R4BP cannot be derived, classical perturbation methods are available to construct the analytical approximate solutions of the periodic motions.The analytical approximate solutions propose a way to inspect the trajectory and exhibit possible parameter design.The most representative work was successfully demonstrated by Richardson [3] who derived the thirdorder analytical solution of halo orbits with small parametric expansions in CR3BP with Lindstedt-Poincaré (L-P) method.The third-order analytical solutions of halo orbits are extensively used as initial conditions in mission trajectories designing progress for numerical integration [4].Similarly, Gómez and Marcote [5] solved a reduced case of the CR3BP, namely, the Clohessy-Wiltshire equation with L-P method.Jorba and Masdemont [6] expanded the Lissajous and halo orbits of CR3BP as formal series of two amplitudes and semianalytically constructed high-order solutions with the assistance of computer.
Similar to CR3BP, it has been analytically proved that the R4BP has some relative equilibrium solutions like the RTBP, and there are families of periodic orbits in the vicinity of equilibrium points [7].The quasi-halo orbits in the Sun-Earth-Moon system were computed from the bicircular model in R4BP by Andreu [8].Papadakis [9] computed the invariant stable and unstable manifolds in R4BP.For the Sun-Jupiter-Trojan asteroid-spacecraft system, the evolution of the families of periodic orbit and their stability around the asteroid and/or the Jupiter in R4BP were calculated by Baltagiannis and Papadakis [7,10].The periodic orbits in the R4BP with two equal masses were explored by Burgos-García and Delgado [11].Recently, Lei and Xu [12] investigated the motions around the equilibrium points of R4BP, in which the motions are expressed as formal series of long, short, and vertical periodic amplitudes.By means of L-P method, series solutions are constructed up to a certain order.
In this study, we model the CR3BP and bicircular case of R4BP by the conservative autonomous system and nonautonomous parametric excited system, respectively, from the vibration point of view, which may shed a new light on the study of periodic orbit investigations.
The classical celestial perturbation methods can give an approximate analytical solution by correcting the amplitude and the period of the linear solutions with the nonlinear terms in the three-dimensional space.However, less attentions have been paid to the relationship of the motions in each dimension of the major and minor axes and their contributions to the nonlinear dynamics.The nonlinear complex modes analysis based on the invariant manifold method by Shaw and Pierre [13][14][15] and Avramov [16] could provide a chance to focus on the nonlinear relations of the motions in different directions.Inspired by Shaw and Pierre and Avramov, the polynomial expansion method for autonomous system and iterative polynomial expansion method for nonautonomous system are innovatively proposed in this paper.By introducing polynomial expansion relations, we transfer the three-degree-of-freedom periodic motions for the nonlinear governing equations into one-degree-of-freedom system for the autonomous CR3BP system.By generating the pseudoautonomous equation for the R4BP system, periodic motions can be obtained by iterative polynomial expansion method.The remainder of this paper is structured as follows.Section 2 introduces the dynamical model R4BP and the CR3BP model.Section 3.1 discusses how to construct high-order series expansions of periodic solution in the CR3BP with polynomial expansion method.In Section 3.2, the iterative polynomial expansion method by transferring the parameter resonance nonautonomous three-degree-of-freedom equations into a pseudoautonomous system is introduced.At last, the conclusions are drawn in Section 4.

Equations of Motion
2.1.Restricted Four-Body Problem.The geometry of the R4BP is shown in Figure 1.The finite masses  1 ,  2 , and  4 and the infinite mass body  3 are denoted as  1 ,  2 ,  4 , and  3 , respectively.The motion of  3 is influenced by the attraction of  1 ,  2 , and  4 , but the motions of the finite masses  1 ,  2 , and  4 are assumed not to be affected by  3 .
In the R4BP,  1 and  2 move about their barycenter perturbed by  4 .The masses  1 ,  2 , and  4 all move about the mass center of the entire system.It is assumed that all bodies are modeled as point masses, and the primaries  1 ,  2 , and  4 move in a plane, with the distance from  4 to the  1 - 2 system being much greater than the distance from  1 to  2 .However,  3 is free to move in all three dimensions.
We define the inertial coordinate system, -, with origin at the center of mass of the four-body system, .The coordinate system, -XYZ, rotates at an instantaneous rate of Ω 1 about the mass center of the four-body system, .The -axis points to the barycenter of  1 and  2 , .The synodic coordinate system is -xyz with origin at .The axis points from the barycenter, , to  2 .-xyz moves with respect to the coordinate system, -XYZ, at an instantaneous rate of Ω.  is the angle between the -axis and the axis, counterclockwise measured.Meanwhile, -xyz rotates with respect to the inertial coordinate system, -, at an instantaneous rate of Nondimensionalization is applied to simplify the form of the equation.The unit length is set as the distance between  1 and  2 .The unit time is set as /2, where  is the period of the  1 - 2 system.We introduce  as the mass parameter of the system; then ( The massive  1 , with mass 1 − , is located at (−, 0, 0) in xyz, and  2 , with mass , is located at The motion of  3 is governed by the following equations with dimensionless time and space with respect to the coordinate system, -.
in which R 1 , R 2 , and R 4 are the position vectors from  1 ,  2 , and  4 to  3 , respectively.And  is the position vector from the system barycenter, , to  3 , with geometrical relation as Note that R  and r are presented in -XYZ and -xyz, respectively, and can be expressed as Apparently, we have time derivative of  with respect to the inertial frame as in which the superscripts "" and "" mean the derivatives are measured with respect to the rotating coordinates -XYZ and -xyz.Combining ( 2) and ( 6), the motion of  3 in synodic coordinate system -xyz is obtained as where Then, we introduce the transformation matrix between -XYZ and -xyz: R 4 can be expressed as Thus, we obtain the completed equation of motion of R4BP in -xyz as where U is the pseudopotential function of the four-body problem and is expressed as

Bicircular Model in Restricted
Four-Body Problem.The bicircular model is a special case of the R4BP and built on the hypothesis that the motions of three primary bodies are restricted.Specifically,  1 and  2 are revolving in circular orbits around their barycenter and the  1 - 2 barycenter, , is moving in a circular orbit around the center of mass of the system .Besides that, the two orbital planes are assumed to coincide with each other [17].
With the aforementioned hypothesis, we consider International Journal of Aerospace Engineering Equation ( 11) is simplified as With expansion of 1/ 4 , equation of motion ( 14) can be further reduced to Since  4 does not affect the relative orbital positions of the other two finite masses, its motion is subjected to Kepler's third law, shown as Therefore, ( 16) equals It is clear that  4 's gravitational force consists of two parts: the direct constant gravity part due the average gravitational force and the parametric excitation part due to the periodic variation of the position.

Circular Restricted Three-Body
Problem.The circular restricted three-body problem (CR3BP) model is a simplification of the R4BP without considering the fourth body's influence,  4 = 0.  1 and  2 still revolve in circular orbits around their barycenter under the influence of their mutual gravitational attractions. 3 moves in the plane defined by the two revolving primaries.The motion of the particle is governed by the following equations with dimensionless time and space units [18]: where Ω is the pseudopotential function of the three-body problem and is expressed as Apparently, all the governing equations of the motions around the libration points in R4BP model, bicircular model, and the CR3BP model belong to the class of gyroscopic systems [19], which always behave as the massless body reaches an extremum in one coordinate as it passes through zero in the other coordinate, such that they are /2 out of phase.

Approximate Analytical Methodology
The governing equations of R4BP and CR3BP in (11) and (19) are widely used in the applications of orbital design, mostly solved by the traditional perturbation methods.In this section, we propose an approximate analytical methodology to derive the polynomial series relations among different directions of the three coupled equations, with which the threedegree-of-freedom problem may be transferred to simple system with one degree-of-freedom.The proposed polynomial relation method can deal with the equation of motions for both autonomous system and the nonautonomous system, which may describe the overall characteristics of the whole system with general rules.

Polynomial Expansion Method for Autonomous System.
First, we describe the polynomial expansion methodology in the CR3BP case.The three second-order nonlinear ordinary differential equations in (19) account for the nonlinear vibrations of the three directions of the conservative system.Since the three equations are coupled by both the gyroscopic terms and nonlinear potential terms, there must exist relations among the three directions.If we can locate the relations, then the three coupled equations can be replaced by one equivalent equation.
Before we implement the routine to find the relations, we need to transfer the reference frame to the principle coordinate system with the origin exact on the libration point, by which the (0, 0, 0) is the equilibrium position and the linear potential terms are decoupled.Now, we transfer the three second-order equations on the principle coordinate system into six first-order equations as where  1 ,  2 , and  3 denote the displacements on the , , and  directions, respectively, and  1 ,  2 , and  3 denote the corresponding velocities.According to the invariant manifold method proposed by Shaw and Pierre [13,20,21], during periodic motion, the positions and velocities of all the directions are functionally related to a single position-velocity pair of the base direction, which we choose arbitrarily, for example, the displacement and velocity on the -axis,  1 and  1 .The position and velocity along the -axis can be assumed as  =  1 and V =  1 , and the positions and velocities on the and -axis are related to  and V. Since we are treating the nonlinear differential equations, nonlinear relations of   and   should be considered as where   and   are real functions of  and V.   and   are the coefficients to be determined.For the weak nonlinear system, the Taylor series expansion of   and   is used to obtain the relations by determining the coefficients of the power series.
In order to determine relations (22), the time derivatives of   and   need to be expanded using the chain rule, as By comparing ( 21) to (23), we can obtain the following equations: where any occurrence of   or   in   is systematically replaced by   or   , respectively.To be detailed, for the libration cases, substituting the Taylor series expansions of   and   into (24), the following equations are obtained: International Journal of Aerospace Engineering Gathering the like terms of the coordinates and velocities provides algebraic equations from which the Taylor series coefficients of   and   can be determined.Finally, substituting the analytical relations of (22) to the -direction equation in (21), we transfer the three-degree-of-freedom problem into one-degree-of-freedom problem, which is a function of .The resulting nonlinear one-degree-of-freedom function of  and its time derivation V can be treated easily by any perturbation methods or numerical methods.The remaining two directions of motion can be expressed as the nonlinear relations of  and V.In the potential applications, the resulting analytical expressions may be as both initial points of computations and constraint conditions for the numerical orbits design.

Iterative Polynomial Expansion Method for Nonautonomous System
. By inspection of (11) or (18), we can concluded that the bicircular model in R4BP is nonautonomous system in the sense of parametric resonance.From the view point of perturbation, when the parametric period is a half of the period of the corresponding autonomous system, the parametric resonance may occur.Since the bicircular model is nonconservative and nonautonomous, the methodology used in the last subsection can not be used directly to the system with time-varying parameters.However, we can substitute the time-varying terms by the base displacementvelocity pair ( 1 ,  1 ), from the mathematical point of view.The iterative idea by Avramov [16] should be considered to fulfill this methodology.
Similarly, we transfer the three second-order equations on the principle coordinate system into six first-order equations. 1 ,  2 , and  3 denote the displacements on the , , and  directions, respectively, and  1 ,  2 , and  3 denote the corresponding velocities.The first step is to derive the nonlinear relations of the and -axis with -axis by determining the coefficients of the Taylor series of (22) by neglecting the timevarying parameters in the governing equations.This step can be carried out as exactly described in the last subsection.However, here the coefficients we derived for the first step are not the final values: we need to optimize these values by iterations.
With the relations based on the autonomous system, which can be expressed as we can transfer the three-dimensional system of  (18) into one-dimensional system on  direction, where higher-order terms other than two have been neglected.The second step is to transfer the time-varying terms as expressions of ( 1 ,  1 ).To study the case of parametric resonance, we only consider the period of the system is approximately a half of the parametric pulsating period, θ ≈ 2, or θ ≈ .
(34) Now, we express the time-varying terms by the functions of ( 1 ,  1 ).Replacing all the time-varying terms by the functions  and  in the governing equations, we can transfer the parameter resonance nonautonomous three-degree-offreedom equations into a pseudoautonomous system.
The last step is to apply the method described in the last subsection again to this set of pseudoautonomous equations by determining the relations of the and -axis with the -axis and transfer the three-degree-of-freedom system into that of one degree-of-freedom.The three-step procedure can be iterated further until satisfactory results are obtained.

Conclusions
In this study, the restricted three-body problem (R3BP) and restricted four-body problem (R4BP) are classified as conservative autonomous system and nonconservative nonautonomous parametric excited system, respectively.From the vibrational point of view, the methodology of polynomial series is proposed to solve for these problems analytically.By introducing the polynomial series relations among the , , and  directions, the three-degree-of-freedom coupled equations are transferred into one degree-of-freedom containing the full dynamics of the original autonomous system for the R3BP.On the other hand, the methodology of polynomial series is combined with the iterative approach to treat the nonautonomous system of bicircular R4BP.The pseudoautonomous equation and the final polynomial series relations and one-degree-of-freedom system can be derived iteratively.
The methodology proposed in this study may be used in the computations of periodic orbits and further be applied in the numerical orbit design as initial conditions and constraints.