#### Abstract

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.

#### 1. 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-moon-moon system, such as the Saturn-Tethys-Telesto-spacecraft, Saturn-Tethys-Calypso-spacecraft, and Saturn-Dione-Helen-spacecraft 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 third-order 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–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.

#### 2. Equations of Motion

##### 2.1. Restricted Four-Body Problem

The geometry of the R4BP is shown in Figure 1. The finite masses , , and and the infinite mass body are denoted as , , , and , respectively. The motion of is influenced by the attraction of , , and , but the motions of the finite masses , , and are assumed not to be affected by .

In the R4BP, and move about their barycenter perturbed by . The masses , , and all move about the mass center of the entire system. It is assumed that all bodies are modeled as point masses, and the primaries , , and move in a plane, with the distance from to the - system being much greater than the distance from to . However, 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, -, rotates at an instantaneous rate of about the mass center of the four-body system, . The -axis points to the barycenter of and , . The synodic coordinate system is - with origin at . The -axis points from the barycenter, , to . - moves with respect to the coordinate system, -, at an instantaneous rate of . is the angle between the -axis and the -axis, counterclockwise measured. Meanwhile, - 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 and . The unit time is set as , where is the period of the - system. We introduce as the mass parameter of the system; thenThe massive , with mass , is located at in -, and , with mass , is located at in -.

The motion of is governed by the following equations with dimensionless time and space with respect to the coordinate system, .in which , , and are the position vectors from , , and to , respectively. And is the position vector from the system barycenter, , to , with geometrical relation as Note that and are presented in - and -, respectively, and can be expressed as

Apparently, we have time derivative of with respect to the inertial frame asin which the superscripts “” and “” mean the derivatives are measured with respect to the rotating coordinates - and -. Combining (2) and (6), the motion of in synodic coordinate system - is obtained aswhere

Then, we introduce the transformation matrix between - and -: can be expressed asThus, we obtain the completed equation of motion of R4BP in - aswhere is the pseudopotential function of the four-body problem and is expressed as

##### 2.2. 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, and are revolving in circular orbits around their barycenter and the - 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 considerEquation (11) is simplified asWith expansion of ,equation of motion (14) can be further reduced toSince does not affect the relative orbital positions of the other two finite masses, its motion is subjected to Kepler’s third law, shown asTherefore, (16) equals

It is clear that ’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.

##### 2.3. 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, . and still revolve in circular orbits around their barycenter under the influence of their mutual gravitational attractions. 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 out of phase.

#### 3. 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 three-degree-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.

##### 3.1. 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 aswhere , , and denote the displacements on the , , and directions, respectively, and , , and 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, and . The position and velocity along the -axis can be assumed as and , and the positions and velocities on the - and -axis are related to and . Since we are treating the nonlinear differential equations, nonlinear relations of and should be considered aswhere and are real functions of and . 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:

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 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 . In the potential applications, the resulting analytical expressions may be as both initial points of computations and constraint conditions for the numerical orbits design.

##### 3.2. 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 displacement-velocity pair (), 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. , , and denote the displacements on the , , and directions, respectively, and , , and 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 time-varying 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 aswe 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 (). To study the case of parametric resonance, we only consider the period of the system is approximately a half of the parametric pulsating period,

Hence, the solution to (27) can be presented in the following form:Substituting (29) back into (27) and using the harmonic balance method, we can obtain a set of five nonlinear algebraic equations with respect to the six parameters (, , , , , and ). By introducing the notationsthe following equations are obtained from (29):The solution of (31) with respect to and is presented in the form of polynomial series:where , are unknown coefficients. Substituting (31) into (32) and collecting the same terms of and letting the coefficients equal zero yield the following linear algebraic equations:where , , , , and is the matrix, the elements of which are expressed in the six parameters .

Solving the two systems of linear algebraic equation (33), polynomial series (32) can be calculated and the following two equations are obtained from these series:

Now, we express the time-varying terms by the functions of (). Replacing all the time-varying terms by the functions and in the governing equations, we can transfer the parameter resonance nonautonomous three-degree-of-freedom 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.

#### 4. 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.

#### Competing Interests

The authors declare that they have no competing interests.

#### Acknowledgments

The authors gratefully acknowledge the support of the National Natural Science Foundation of China (NNSFC) through Grants nos. 11402007, 11290152, 11672007, and 11322214 and the Funding Project for Academic Human Resources Development in Institutions of Higher Learning under the Jurisdiction of Beijing Municipality (PHRIHLB).