Comparisons between the circular restricted three-body and bi-circular four body problems for transfers between the two smaller primaries

Important properties of the dynamics of a spacecraft can be obtained from the Circular Restricted Three Body Problem and the Bi-Circular Bi-planar Four Body Problem. In this work, both systems are compared under the perspective of the costs involved in a transfer between the smaller primaries. An analytical approach shows several properties of the perturbation due to the gravity of the Sun and the motion of the smaller primaries around it over a spacecraft in the region of interest, like its behavior at and around the barycenter or at any point in a circle around the Sun. The costs involved in transfers between the smaller primaries are numerically evaluated and analyzed using the newly developed Theory of Functional Connections. The results show that the influence of this perturbation over the costs is significant for systems like the Sun–Earth–Moon or Sun–Mars–Phobos. On the other hand, it is also shown that this influence may be negligible for other very different systems, like the Sun–Saturn–Titan or Sun–Ida–Dactyl. Maps of perturbation are drawn in the region of interest, which can be used for mission designers. Finally, a new approach to describe the influence of the Sun over the tides of the smaller primaries is proposed under the Four Body Problem model.

The three-body problem is described by Newton's equations of motion for three bodies subjected to their mutual gravitational attraction. The circular restricted three-body problem (CR3BP) is a particular case where one of the bodies has a negligible mass (e.g. a satellite) and the motion of the other bodies are circular around their barycenter. Although there is no analytical closed (other than series) solution to this problem 1 , this model has been proven to be very useful in astrodynamics to find general features of orbits. This is possible due to its simplification from many perturbations to the few main ones in the region under consideration.
A combination of two constrained "three-body problems" was structured by 2 to find qualitative properties of the motion of satellites under the gravitational influence of the Earth, Moon, and Sun. In general, the massive bodies in the first restricted three body problem are a planet (or an asteroid) and its moon, while the massive bodies in the second restricted three body problem are given by the Sun and the other two pairs combined. In a subsequent paper, mathematical descriptions of the very restricted four body problem was done by the same author in 3 . In comparison with the CR3BP, this model takes into consideration the gravitational influence of the Sun by adding a term which is a function of the position of the Sun relative to the Earth-Moon system.
The existence of periodic motion around the classical libration points of the CR3BP was shown for this four body problem (4BP) in 4 . Since then, the four body problem has been successfully investigated for specific masses and positions of the bodies in the space (under symmetrical configurations) [5][6][7][8] . Nevertheless, the restricted four body problem has attracted the interest of many researchers in astrodynamics, because it can be used to find properties of the motion of a satellite-its mass is much lower than the other three bodies. The restricted four body problem has also been investigated under specific configurations as the one in which all the primaries are located in a straight-line equilibrium configuration 9 or all the primaries are located in the vertices of an equilateral triangle with equal 10 or different masses 11 -a problem similar to the Hill's approximation in which the mass of one of the primaries is much lower than the mass of the others. Furthermore, the conditions in which

Mathematical models
In this section, the mathematical tools and the two models adopted in this research are shown and compared.
The bi-circular bi-planar four body problem. A vector h * connects the center of an inertial frame of reference to the center of a rotating frame of reference. The Sun is located at the center of the inertial frame, and the barycenter of a M e -moon system is located at the center of the rotating frame of reference. The inertial frame of reference is denoted by a star (*), for clarity purpose. Thus, the position of a particle in the rotating frame of reference ( r ) with respect to the inertial one ( r * ) is The acceleration can be written as 28 where ω is the angular velocity of the rotating frame. The derivatives with no star ( d dt ) are taken in the rotating frame, i.e. it does not take into consideration the motion of its base, while the derivatives with star ( d * dt ) are taken in the inertial frame. Hence, the equation of motion of a particle written in the rotating frame of reference is where a is the acceleration acting over the particle.
In the bi-planar, bi-circular four-body problem, the body M e and the moon rotate in a circular orbit around their barycenter, which in turn rotates in a circular orbit around the Sun. Moreover, the two orbital planes coincide with each other. The two frames of reference can be seen in Fig. 1.
The barycenter of the M e -moon rotates in a circular motion around the Sun, according to where R s is the constant distance between the Sun and the barycenter, and θ 1 = (ω r t + γ − π) is the angle between h * and the positive side of the x * axis, where γ − π is an initial phase and ω r is the angular velocity of h * . The second derivative of h * with respect to the inertial frame, written in the rotating frame of reference is where µ s is the gravitational parameter of the Sun and ω s = ω r − ω.
The vectors that locate the spacecraft with respect to M e , the moon, and the Sun are r e , r m , and r s , respectively. The first two are given by r e = r + d ex and r m = r − d mx , where x is an unitary vector along the x axis, d e = Rµ m /(µ e + µ m ) and d m = Rµ e /(µ e + µ m ) are the distances between the barycenter and M e and the barycenter and the moon, respectively, where µ e and µ m are the gravitational parameters of M e and the moon, respectively. The last one is given by where R s = R s cos θ 2 , sin θ 2 , 0 locates the Sun in the rotating frame, where θ 2 = θ 2 (t) = ω s t + γ.
Hence, the equation of motion (3), according to the bi-circular bi-planar four body problem, becomes The circular restricted three body problem. In the case where the Sun is not taken into consideration and the center of the inertial frame of reference is located at the barycenter of the M e -moon system, the equation of motion becomes The perturbation due to the gravity of the moon is It should be noted that the motion is Keplerian when this term is null in Eq. (8).
(4) h * = R s cos θ 1 , R s sin θ 1 , 0 cos θ 2 , sin θ 2 , 0 .  www.nature.com/scientificreports/ Perturbation due to the gravity of the Sun. Note that, in the case where µ s = 0 , Eq. (7) is reduced to the planar circular restricted three-body problem for the M e -moon-satellite system given by Eq. (8), which does not take into consideration the influence of the Sun. Thus, the perturbation of the Sun over the satellite is given by the last two terms of Eq. (7), which can be written as Note that this perturbation does not represent only the acceleration due to the gravity of the Sun, it is a pseudo specific force instead. The second term on the right side comes from the fact that the rotating frame of reference (and also the pair M e -moon) is in an accelerated motion around the Sun. On the other side, it can be seen as a perturbation in the CR3BP due to the presence of the Sun.
Orbit transfer using a bi-impulsive maneuver. Comparisons between both models will be shown for several bodies in the Solar system, with different masses and distances. This is done through an important subject of the astrodynamics, which is the costs associated to a transfer of a spacecraft from an orbit around a massive body to another orbit close to its moon. The transfer is done in the following way: a spacecraft is initially in a circular orbit of radius r 0 around the massive body M e , when a first impulse is applied to it at a point A of this orbit. This impulse is applied such that the spacecraft travels to a second point B close to the moon in a given time of flight. The distance between B and the moon is ρ 0 . A second impulse is applied at the point B such that the spacecraft enters in a circular orbit around the moon. The magnitude of the differences of the velocity of the spacecraft before and after the first impulse at point A plus the magnitude of the difference of the velocity of the spacecraft before and after the second impulse at point B is the cost of the transfer.
This problem is known as the "two point boundary value problem", i.e. the maneuver is constrained such that the position of the spacecraft at the initial time is A and its position at the final time is B. The difference between the final and initial times is the time of flight, which is assumed to be a given parameter of the transfer. The solution of this problem is obtained using the Theory of Functional Connections 29,30 . This is a very efficient approach to solve both linear and and nonlinear differential equations 31,32 subjected to two constraints in the positions-the "two point boundary value problem". More details in the use of the Theory of Functional Connections to solve the orbit transfer problem can be seen in 20 . The numerical evaluations was done using the Python language combined with a TFC Python module freely available in 33 . An algorithm is used to find local minima around the neighborhood of a point in the space of parameters, which are varied each at a time. The costs in this neighborhood are compared and the parameters associated to the lowest value of the cost are selected as the optimal local maneuver. This is done for all of the parameters, several times. The step is reduced for each successive iteration, until the variation of the cost after each step is no more than 10 −3 m/s , which is assumed to be our acceptable error. In this way, the minimum cost associated to a given time of flight is found through the algorithm described above.
The cost associated to the time of flight is the minimum for the given transfer time, and it is found for each of the two models: the CR3BP and the bi-circular bi-planar 4BP. In each case, the data are fitted using Legendre's polynomials with its 20 first terms. The resulting curves are the costs as function of the time of flight for the 3BP (named cost 3BP ) and the 4BP (named cost 4BP ). The cost gain is defined as the cost required by the 3BP minus the cost required by the 4BP, given by Thus, the cost gain shows the gain of the 4BP with respect to the 3BP in terms of variation of velocity required by the maneuver. The relative cost gain is defined as the cost gain divided by the cost associated to the 3BP, according to Thus, this ratio compares the cost gain with the cost of the transfer. Finally, the relative cost gain per time of flight is defined as the relative cost gain divided by the respective time of flight: This is an important factor, because longer times of flight tend to accumulate the effect of the perturbation over the costs, and this index tends to cancel this dependency.

Results
The results obtained from simulations using the two models described above are shown next.
Perturbation at the barycenter. In the case where the spacecraft is located at the barycenter, its position is given by r = 0 , where 0 is the null vector. Hence, according to Eq. (6), its position with respect to the Sun is r s = −R s . Using this result, the perturbation given by Eq. (10) becomes p s = 0 . Hence, there is no perturbation from the Sun at the barycenter of the Earth and the Moon-both models are identical at this point.
Perturbation on a sphere around the Sun. The perturbation given by Eq. (10) can be written as In this case, the perturbation due to the Sun is directed toward the barycenter and its magnitude is proportional to the displacement.
Perturbation around the barycenter. The magnitude of the first term on the right side of Eq. (14) is much smaller than each of the last two terms. On the other side, the last two terms tend to cancel each other for r s close to R s , which is the case considered now. Hence, neither term can be neglected. A series expansion of Eq. (10) around the position of the barycenter 0, 0, 0 -neglecting terms of order greater than one and crossed terms on x, y, z-is given by Hence, its magnitude is For a constant z, the dominant term inside the square root is (5x 2 + 5y 2 ) . Thus, constant values of the magnitude of the perturbation are distorted concentric circles around the barycenter. This distortion depends on the value of θ 2 , or the time, because θ 2 is a function of time. In the space, the term 2z 2 is taken into account, and constant values of this magnitude are distorted spherical shells around the barycenter, and it tends to increase faster with respect to x and y in comparison to its dependence with z.
Relative perturbation. The relative perturbation is defined as the ratio between the magnitude of the perturbation and the magnitude of the gravitational attraction of the moon, thus it is given by p s / p m . This ratio shows a very important relationship between these two variables, because p m is the perturbation due to the moon over a Keplerian orbit, and p s is the perturbation over the orbit obtained using the CR3BP. The larger this number, the more important is the perturbation due to the Sun in comparison with the perturbation due to the moon over a Keplerian orbit. The relative perturbation is negligible around the barycenter (because �p s � = 0 at the barycenter) and close to the moon (because p m is much larger when the position tends to r ≈ d mx ). On the other hand, this quantity depends on several parameters between the barycenter and the moon.
The Sun-Earth-Moon system. The results are shown below in this section for the values of the parameters of the Sun-Earth-Moon system shown in Table 1.
The direction of the perturbation p s is represented by the streamlines shown in Fig. 2, as function of the coordinates, for , and θ 2 = 270 • (down right). The magnitude of the perturbation p s is shown by the gradient color, from dark purple for lower values of the magnitude of the perturbation to brighter yellow for larger values. Although the angle θ 2 is a function of time ( θ 2 = θ 2 (t) = ω s t + γ ), the effects on the perturbation of its oscillation from 0 to 2 π can be seen by the drawn of the discrete sequence of its values shown in Fig. 2. An analogous drawn is done in Fig. 3 for the perturbation, as a function of x and z, for y = 0 . The direction of the Sun is perpendicular to the x-z plane in the case where the relative position of the Sun is θ 2 = 90 • or θ 2 = 270 • . In these cases, the direction of the perturbation is toward the barycenter. Note that, in the regions where the position is such that r s ≈ R s , the perturbation is approximated by the same value as the one it has on a sphere around the Sun, shown in the subsection above, i.e. p s ≈ −µ s r/R 3 s . This behavior can be seen in Fig. 2 for the regions close to a straight line that crosses the barycenter and is perpendicular to the direction of the Sun. Furthermore, this behavior can also be seen in Fig. 3, in the case where θ 2 = 90 • or θ 2 = 270 • . Note that the sphere given by the (14) 3 cos(2θ 2 ) x 2 − y 2 + 6xy sin(2θ 2 ) + 5x 2 + 5y 2 + 2z 2      Table 2. The relative perturbation in the plane x − y can be seen in Fig. 6. Note that, in this case, the barycenter of the pair Mars-Phobos is almost in the center of Mars, far from its surface. On the other hand, the magnitude of the perturbation due to the gravity of Phobos is very small. Hence, although the distance Mars-Phobos is small and the region around the pair is close to the barycenter (where p s is null), the relative perturbation is one order of magnitude higher in comparison with the Earth-Moon case. Thus, the motion described by the CR3BP converges with the one described by the bi-planar bi-circular 4BP only in a small region around Phobos. This comparison also helps to understand the importance of each body in the dynamics. The minimum cost as a function of the time of flight is shown in Fig. 7 (up left part) for a transfer from an initial circular orbit around Mars of 167 km of altitude (above the Mars radius of 3389.5 km) to a final circular orbit around Phobos with 10 km in altitude (21.2667 km from its center). The cost gain is shown in the up right    The Sun-Saturn-Titan system. In this subsection, the relative perturbation is shown in Fig. 8 for the Sun-Saturn-Titan system, whose parameters can be seen in Table 3. The relative perturbation is also presented for several values of the relative position of the Sun given by θ 2 . Note that it is about one order of magnitude lower than the value for the case of the Earth-Moon system and almost two orders of magnitude lower than the Mars-Phobos case.
The cost as a function of the time of flight is shown in Fig. 9-up left side-from a circular orbit with semimajor axis of 100 × 10 6 m to another circular orbit around Titan with 100 km of altitude (semi-major axis equals to 2.657473 × 10 6 m ). The cost gain, the relative cost gain, and the relative cost gain per time of flight are also shown in this figure. Note that the distance Saturn-Titan is much larger than the distance Earth-Moon. Thus, the region of interest (which is the region where the transfer is realized) between these two main bodies is further from the barycenter, where the magnitude of the perturbation ( p s ) is null. On the other side, Titan has more mass than our Moon, which makes its cost contribution to the transfer more significant. Thus, the relative perturbation for this case is lower than in the Earth-Moon system, in the respective regions of interest. The lower order of magnitude of the average relative index cost gain per time of flight can explain the lower relative perturbation of the Saturn-Titan in comparison with the Earth-Moon system (compare Fig. 4 with 8 and Fig. 5 with 9).
The Sun-Ida-Dactyl system. The relative perturbation is shown in Fig. 10 for the Sun-Ida-Dactyl system, whose values of the parameters are shown in Table 4. The relative perturbation is two orders of magnitude lower than the similar values in the Earth-Moon system in the region of interest around the bodies, where the path of a short-time transfer is located. Note that the distance Ida-Dactyl is only 90.5 km, according to Table 2. The magnitude of the perturbation ( p s ) is very low for this system, due to the proximity of this region to its barycenter (close to the center of Ida), which can explain the very low relative perturbation in the region around the bodies, about two orders of magnitude lower than those for the Earth-Moon system evaluated in the respective regions.  Fig. 11, for a transfer from a circular orbit around Ida whose radius is 60 km, to another circular orbit around Dactyl whose radius is 2 km. Clearly, all indices show that both the 3BP and the 4BP coincide for the Ida-Dactyl system. As explained before, the closer the region is to the barycenter, the more the results for the 3BP coincide with the ones obtained from the 4BP. The average value of the relative cost gain per time of flight is only 0.001% per day, which is the lowest one among the cases shown in this section. Note that the numerical proceedings was done with higher precision for this case, in order to evaluate the very small differences (the average relative cost gain is www.nature.com/scientificreports/ only 0.0005%). Possibly, a more refined numerical evaluation could give an even smaller difference between the 3BP and the 4BP, at the cost of an increased computational cost. On the other side, the smallest relative cost gain per time of flight is in agreement with the smallest relative perturbation among the cases studied in this section.
Comparisons of the results obtained for the Solar system. In this section, the results obtained in the previous sections for several pairs in the Solar system will be analyzed, based in the data obtained from the previous subsections shown in Table 5, where the averages of the indices defined in "Orbit transfer using a biimpulsive maneuver" section are shown. The values of the relative cost gain per time of flight may be separated into two groups, whose values are of the same order of magnitude. The first one is given by the Sun-Earth-Moon and the Sun-Mars-Phobos systems and the second group is given by the Sun-Saturn-Titan and Sun-Ida-Dactyl systems. Although the gravitational parameters of Phobos is much smaller than the gravitational parameter Table 3. Values of the parameters for the Sun-Saturn-Titan system 35 . This means that the spacecraft is much closer to the barycenter of the smaller primaries in a Mars to Phobos transfer than it is in the case of an Earth to Moon transfer. The final result is that the savings are similar for transfers in both systems. The comparison between the Sun-Saturn-Titan and the Sun-Ida-Dactyl can also be done using a similar analysis. On one side, the gravitational paramater of Titan is ten orders of magnitude greater than the gravitational parameter of Dactyl, and, on the other side, the parameter R is five orders of magnitude lower www.nature.com/scientificreports/ for the Sun-Ida-Dactyl in comparison with the Sun-Saturn-Titan system. This means that travels in the Ida-Dactyl system are made much closer to the barycenter compared to travels between Saturn and Titan. The final effect is that, for both systems, the influence of the perturbation due to the Sun over the fuel savings is negligible.
Variation of the parameters for the Solar system. In this section, the behavior of the relative perturbation p s / p m is analyzed as a function of the parameters involved in the problem. Note that the relative perturbation is a function of x, y, z, and θ 2 . In order to simplify the analysis, two averages are performed. The first average is taken for the relative perturbation over a straight line from the position of µ e (−d e , 0, 0) to the position of µ m (d m , 0, 0) , according to  Once these averages are taken, the dependency of the relative perturbation on the position on the space (x, y, z) and on the direction of the Sun (θ 2 ) are removed. The average on the position provides a good estimation of the order of magnitude of the relative perturbation. It was seen that the magnitude of the perturbation p s is a distorted sphere around the barycenter, thus, the average on the direction of the Sun θ 2 can also provide a good estimation of its order of magnitude. Furthermore, the major interest here is to investigate the general behavior of the relative perturbation, not its local variation. Assuming that µ s is the gravitational parameter of the Sun, the averages are functions of four parameters: µ e , µ m , R s , and R.
The first case to be analyzed is the dependency of the averages on the parameter µ e . Note that the parameter R is, in general, lower than the sphere of influence of µ e . Hence, for this case, the value of the parameter R is R = R SOI /2 , where the sphere of influence is R SOI = R s µ e µ s 2 5 , according to 39 . In this case, the gravitational parameter of the moon is µ m = µ e × 10 −5 . Finally, the distance from the Sun is R s = 1 AU . In this case, the averages as functions of µ e are shown in Fig. 12.
The second case to be analized here is the dependency of the averages on the parameter µ m . The other parameters are such that R s = 1 AU , R = R SOI /2 , and µ e is the gravitational parameter of the Earth, whose value is shown in Table 1. In this case, the averages as functions of µ m are shown in Fig. 13.
The variations of the averages with the parameter R s are also analyzed. In this case, the other parameters are such that R = R SOI /2 , µ e is the gravitational parameter of the Earth, shown in Table 1, and µ m = µ e × 10 −3 . In this case, the evaluations showed that they have averages �p s � �p m � = 1.42144 for every value of R s in the range 10 9 m < R s < 10 13 m. Table 5. Values of the parameters and the indices for several systems, where AV 1 is the average cost gain, AV 2 is the average relative cost gain, and AV 3 is the average relative cost gain per time of flight.  Fig. 14. The other parameters are given by R s = 1 AU , µ e is the gravitational parameter of the Earth, shown in Table 1, and µ m = µ e × 10 −3 .

Conclusion
The present paper developed a new method to accurately measure the differences between the CR3BP and the bi-circular bi-planar 4BP dynamical models, which is particularly important in the regions of transfers around the smaller primaries. This method can be used to choose the model to be used in real mission design, evaluations of trajectories, optimization of computational time and fuel consumption, with high accuracy.
To develop this method, a term provided by the bi-circular bi-planar four body problem due to the influence of the Sun is added. It is seen as a perturbation, which explicit shows the differences between the bi-circular bi-planar four body problem and the CR3BP. The lower this term, the closer the bi-circular bi-planar four body problem is to the CR3BP.
Firstly, it was analytically shown (for any mass of each one of the three bodies) that the perturbation of the Sun is null at the barycenter of the system composed by M e and its moon. After that, it was also shown that the perturbation is linearly proportional to the position of the satellite with respect to the barycenter of the M e -moon system, when the satellite is at a fixed distance R s from the Sun, which can be seen as a straight line transversal to the direction of the Sun, for a satellite close to the Earth. Finally, the numerical results confirmed the concentric oval curves for the same values of the magnitude of the perturbation, although it also shows that the direction of the perturbation have different symmetries.
In general, the results show that the closer the spacecraft is to the Earth-Moon barycenter, the lower the magnitude of the perturbation. It happens due to the inclusion of the Sun in the equations of motion, and, hence, the closer the bi-circular bi-planar four body problem is to the circular restricted three body problem.
Investigations were performed by evaluating orbit transfers in several pairs of bodies in the Solar system, like Earth-Moon, Mars-Phobos, Saturn-Titan, and the Ida-Dactyl systems. The results showed that the lowest perturbation of the Sun is for the Ida-Dactyl pair, which means that the CRTBP and the bi-planar bi-circular  www.nature.com/scientificreports/ 4BP gives very similar results for this system. The system Saturn-Titan showed to have the second lowest differences between the two models. The third lowest differences were found for the Earth-Moon system. Finally, the larger differences were found for the Mars-Phobos system. It means that the CR3BP is not a good model to study this system. The influences of these costs are strictly linked to the fuel savings when using the bi-circular 4BP in comparison with the CR3BP, depending on the position of the Sun relative to the smaller primaries. These savings are evaluated and shown in this paper considering short time transfers. Although the coincidence of the CR3BP and the bi-planar bi-circular 4BP depends on the combinations of the masses of the pair M e and its moon, its internal distance, and their distances from the Sun, it was found that the differences between the models increase with the masses of the pair and with their internal distance from the common barycenter.