Periodic Orbits of the Restricted Three-Body Problem Based on the Mass Distribution of Saturn’s Regular Moons

: This paper uses the Kolmogorov–Smirnov test to perform a ﬁtting analysis on the mass data of Saturn’s regular moons and found that the lognormal distribution is its best-ﬁtting distribution with an extremely high p -value of 0.9889. Moreover, novel dynamic equations for the variable-mass restricted three-body problem are established based on the newly discovered distribution of mass data, rather than the empirical Jeans’ law, and the Lindstedt–Poincaré perturbation method was used to give the approximate analytical periodic orbits near the Lagrangian point L 3 . Furthermore, this paper also discusses the inﬂuence of the three-body gravitational interaction parameter, the variable-mass parameter of the third body, and the scale parameter in the statistical results on the periodic orbits and the position of the Lagrangian point L 3 through numerical


Introduction
Humanity will be challenged by the threat of resource shortages in the future. Some small celestial bodies, planets, and moons may become a breakthrough for obtaining resources. In October 2019, the Sheppard team of the Carnegie Institution in the United States announced the discovery of 20 new moons of Saturn [1]. To date, the Saturn family has 82 members, making it the planet with the largest number of moons in the solar system (see [2] for specific classification). These natural satellites are divided into regular and irregular moons based on their distance from Saturn, orbital eccentricity, and inclination. From their dynamic behavior and the composition and structure of matter, the Saturnian system is similar to a miniature solar system. For example, scientists discovered an atmosphere and liquid water on Titan. Some scientists have thought that the moon might be an asteroid captured by its host star or the debris produced by its parent body or other moons colliding together [3]. The sources of collisions are diverse, such as main belt asteroids and various comets [3,4]. Dorofeeva [5] believed that the original matter that composes the regular moon Enceladus came from a particular comet, and Castillo-Rogez et al. [6] considered that Phoebe came from a library of C-type asteroids. Gao et al. [7,8] conducted a fitting distribution study on the irregular moon distribution of Jupiter. According to the p-value obtained by the Kolmogorov-Smirnov (K-S) test, they proposed that the log-logistic distribution is the best distribution of their main physical characteristics. The best distribution model inferred from the data can help predict unknown moon physical properties and even help humans discover new natural satellites.
In celestial mechanics, the restricted three-body model describes the motion of an object with a relatively infinitely small mass under the gravitational force of two finitemass bodies (that is, the main star). It is called the circular restricted three-body problem if the finite mass bodies move in a uniform circular motion around their mass center. This model can approximate the local dynamics of planets or moons with small orbital eccentricities. Considering that the mass and radiation of celestial bodies are constantly

Best-Fitting Distribution of the Regular Moons' Mass
To explore whether the mass data of Saturn's regular moons (see Table 1) have statistical regularity, this section uses the single-sample K-S test to evaluate the mass data of the regular moons. The K-S test is a method to test whether a certain set of sample data comes from a specific distribution. The main idea of this test method is to compare the empirical distribution function of the sample data with the selected theoretical distribution function and check whether the sample data come from this theoretical distribution function according to the difference between the two. The main steps of the test are as follows: first, set the null hypothesis H 0 : the sample data obey a particular distribution, then the alternative hypothesis H 1 : the sample data do not obey the particular distribution. Suppose we use D = max |F n (x) − F 0 (x)| to represent the maximum value of the gap between the empirical distribution function F n (x) and the theoretical distribution function F 0 (x) of the sample data. In that case, the critical value D(n, α) is determined by the number of known sample data n and the significance level α, and if the value of D falls within the rejection range, the null hypothesis is rejected. Otherwise, the null hypothesis is not rejected. In the hypothesis-testing problem, relative to the D-value, the p-value (0 ≤ p ≤ 1) can also reflect the credibility of statistical inference. The p-value is the minimum significance level that can reject the null hypothesis H 0 in the hypothesis-testing problem. In other words, the p-value is the probability of falsely rejecting the null hypothesis H 0 , which is called a Type 1 error. The null hypothesis may become more significant as the p-value increases and decrease as the p-value decreases. When the p-value is smaller than the selected significance level α (this article selected α = 0.05 according to international practice), the null hypothesis H 0 is considered invalid. That is, the sample data do not obey the given theoretical distribution. When multiple distribution functions can fit the same sample data, the distribution function with the most significant p-value can be used as the best-fitting distribution of the sample data. For the mass data of Saturn's regular moons, a total of 22 common distribution functions were fit in this section, and the specific fitting results are shown in Table 2. It is not difficult to find that the maximum p-value can reach 0.9889; the corresponding best-fitting distribution is the lognormal distribution, and its probability density function (PDF) is: the PDF of the normal distribution corresponding to this distribution is: Here, the location parameter −∞ < µ < +∞ is the mean of the logarithmic value, and the scale parameter σ ≥ 0 is the standard deviation of the logarithmic value. The lognormal distribution is sometimes called the Galton distribution, which can be obtained after logarithmic transformation of the normal distribution. It should be noted that the lognormal distribution requires that each datum X be a positive number because ln(X) is meaningful only when X > 0. If X follows the lognormal distribution with parameters µ and σ, then ln(X) follows the normal distribution with mean µ and standard deviation σ. Similarly, if X follows the normal distribution with parameters µ and σ, then exp(X) obeys the lognormal distribution. µ refers to the position parameter of the lognormal distribution, not the position parameter of Saturn's regular moons. From the results of statistical inference, the estimated value of parameter µ in Equations (1) and (2) is 39.1272 with the corresponding confidence interval being [35.7438, 42.5106], and the estimated value of parameter σ is 7.82406 with the corresponding confidence interval being [6.05109, 11.0738]. The cumulative distribution function (CDF) of moon mass data and the CDF corresponding to the best-fitting distribution are shown in Figure 1. It is easy to find that the lognormal distribution can fit the mass distribution of Saturn's regular moons well.

Dynamic Model with Variable Mass
Based on the law of mass distribution in the previous section, this section mainly analyzes the dynamic behavior of the restricted three-body problem composed of the Sun, Saturn, and the third body with variable mass (the mass is negligible relative to the Sun and Saturn). The third body here can be regarded as the "parent planet" of Saturn's moons because it has long been suspected that the moons were formed when the larger bodies disintegrated. According to the results of statistical inference, let the mass m of the third body obey the lognormal distribution with the location parameter µ and the scale parameter σ; then, ln(m) obeys the normal distribution with mean µ and standard deviation σ. Note that if the mass m of the third body is a function that changes with the evolution time τ, then its probability density expression is: Performing the time scale transformation t = (ln m(τ) − µ) 2 /2, we have: then: We considered the movement of the third body in the gravitational field of the Sun and Saturn and established a center-of-mass rotating coordinate system O-UVW: taking the center-of-mass of the Sun-Saturn system as the origin, the plane of relative motion as the orbital plane, the U axis passing through the line of the Sun and Saturn, the V axis perpendicular to the U axis in the plane of motion, and the W axis passing through the center-of-mass of the Sun-Saturn system and rotating perpendicular to its plane of motion. Normalize related physical quantities. Let the mass of the Sun be m 1 and the mass of Saturn be m 2 ; then, the mass ratio parameter µ = m 2 /(m 1 + m 2 ) ≤ 1/2, m 1 = 1 − µ, m 2 = µ. Therefore, the dimensionless dynamic equation of the third body can be described as: where the potential energy function is: where k is the three-body interaction parameter, ρ 2 1 = (u + µ) 2 + v 2 + w 2 and ρ 2 2 = (u − (1 − µ)) 2 + v 2 + w 2 are the distances from the third body to the Sun and Saturn, respectively, and q 1 and q 2 are the radiation factors of the Sun and Saturn, respectively. For more explanation, please refer to [15,23,24]. To simplify Equation (6), according to the space-time transformation formula: where γ = m/m 0 represents the mass ratio and m 0 is the mass of the third body at the initial moment. Substituting Equations (5) and (8) into Equation (6) and letting q = 1/2, l = 0, Equation (6) becomes:ξ where the potential energy function is: here

Periodic Orbits near the Lagrangian Point L 3
Without loss of generality, we analyzed the periodic orbits around the Lagrangian point L 3 . Consider the dynamic equations of the third body on the U-V plane: where: By introducing the transformation ξ = ξ L 3 + x and η = y, where (ξ L 3 , 0) is the coordinate of Lagrangian point L 3 , the dynamic Equation (11) becomes: Using Taylor's formula to expand the right side of Equation (14) to the second-and third-order, respectively, we obtain: The coefficient expressions of the expansion items above are listed as follows: Assume that the solutions of Systems (15) and (16) with respect to the perturbation parameter ε(|ε| ≤ 1) are expressed as: (17) Substituting Equations (17) and (18) into Equations (15) and (16), respectively, and comparing the degree of parameter ε, we have: Note that Hu [34] used the Lindstedt-Poincaré perturbation method to give its periodic solution of Equation (19) in the following form: x 1 = a 1 cos(ωt) + b 1 sin(ωt), y 1 = a 2 cos(ωt) + b 2 sin(ωt), (22) where ω = 2π/T and T is the period. Substituting Equation (22) into Equation (19) yields: (K 1 + ω 2 )a 1 + 2ωb 2 cos(ωt) + (K 1 + ω 2 )b 1 − 2ωa 2 sin(ωt) = 0, If Equation (19) has a nonzero periodic solution, it only needs to satisfy the coefficient determinant: Selecting the parameters a 1 = 1, b 1 = 0, correspondingly a 2 = 0, b 2 = −2ω/(P 1 + ω 2 ), then the periodic solution of Equation (19) is: Similarly, the periodic solution of Equation (20) can be written as follows: x 2 = a 3 + a 4 cos(ωt) + a 5 cos(2ωt) + b 3 sin(ωt) + b 4 sin(2ωt), y 2 = a 6 cos(ωt) + a 7 cos(2ωt) + b 5 sin(ωt) + b 6 sin(2ωt), then it can be obtained by a similar method that: x 2 = a 3 + a 5 cos(2ωt), where: Therefore, the periodic solution of the second-order expansion of the system (15) at the Lagrangian point L 3 with respect to the perturbation parameter ε is: In the same way, the periodic solution of the third-order expansion of system (16) at Lagrangian point L 3 concerning the perturbation parameter ε is: where ω is the positive real root of the equation (24), and:

Numerical Simulation
This section discusses the influence of the scale parameter σ of the regular moons' mass distribution, the three-body interaction parameter k, and the variable-mass parameter γ on the dynamic behavior of the periodic orbit near the Lagrangian point L 3 . According to Stefan-Boltzmann's theorem [35], the radiation factor parameters of the Sun and Saturn are q 1 ≈ q 2 ≈ 1. In addition, the actual mass ratio parameter of the Sun and Saturn is µ = 0.0002857.

Influence of the Scale Parameter σ
According to the value range of parameter k in [23], we chose k = −0.03, and γ = 0.5, ε = 0.1. When the scale parameter σ of the lognormal distribution is 1, 2, 7.82406, and 100, the second-and third-order approximate periodic orbits near the Lagrangian point L 3 on the x-y plane are shown in Figure 2. The degree of agreement between the second-order and third-order approximate periodic orbits increases as σ increases. To visualize the influence of the change in the value of σ on the periodic orbit, we depict the third-order periodic orbits corresponding to different values of σ in the same three-dimensional coordinate system and project them to the three coordinate axis planes, as shown in Figure 3.  When the scale parameter σ located on the denominator of the potential energy function decreases, we found that the shape of the periodic orbits becomes complicated, as shown in Figures 2 and 3. With the increase in σ in Figure 3, the shape of the third-order approximate periodic orbit hardly changes, the position of Lagrangian point L 3 gradually moves away from the center-of-mass and then stays near (−0.7018, 0), and the periodic orbit in the directions of x and y finally stabilizes in the area of (−0.8, −0.6) × (−0.2, 0.2). The range of σ is limited to the confidence interval [6.05109, 11.0738]. In this case, we found that not only the second-and third-order periodic orbits are very close to each other under a given σ-value in the confidence interval, but also the shape and amplitude of the orbit, and the position of the Lagrangian point L 3 remains almost unchanged. This also indicates that when the moons' mass obeys the best-fitting distribution, the statistical error has a minimal effect on the periodic orbit near the Lagrangian point L 3 .

Influence of the Three-Body Interaction Parameter k
We selected the parameters γ = 0.5, σ = 0.65, and ε = 0.1. The approximate periodic orbits near the Lagrangian point L 3 on the x-y plane are shown in Figure 4 when k takes −0.1, 0, 0.3, and 0.7. To examine the influence of the change in k on the periodic orbit and the position of the Lagrangian point L 3 , we describe the third-order periodic orbits corresponding to different values of k in the same three-dimensional coordinate system and projected them to the three coordinate axis planes, as shown in Figure 5.
From Figures 4 and 5, we found that the shapes of the second-and third-order periodic orbits have slight differences when k takes a given value. At this time, the degree of agreement between approximate periodic orbits of different orders cannot be effectively improved by increasing the orbital order. Moreover, as shown in the projections of Figure 5, the periodic orbit amplitudes along the x-direction do not change with the increasing k, but they gradually increase along the y-direction. In addition, the periodic orbits and the position of the Lagrangian point L 3 move away from the center-of-mass. Similarly, to explore the influence of the variable-mass parameter γ on the periodic orbits and the position of the Lagrangian point L 3 , we describe the third-order periodic orbits corresponding to different values of γ in the same three-dimensional coordinate system and projected them onto the three coordinate axis planes, as shown in Figure 7.  Figures 6 and 7 show that the difference between the second-and third-order approximate periodic orbits is considerable when γ is small. This is because a smaller γ corresponds to a smaller third-body mass; in addition to being affected by the gravitational force of the Sun and Saturn, it is also vulnerable to the gravitational or perturbation factors of other celestial bodies, resulting in poor orbital stability, and there is an apparent difference between orbits of different orders. When γ is close to one, the degree of agreement between the second-and third-order periodic orbits is near perfect, and the second-order periodic orbit can be used to approximately replace the third-order periodic orbit. As the mass of the third body decreases, that is as γ decreases, the periodic orbit will shift towards the center-of-mass, and the position of the Lagrangian point L 3 will gradually approach the center-of-mass of the system. However, the change in γ hardly affects the amplitude of the periodic orbit, and the areas of the periodic orbits are also approximately equal.

Conclusions
Based on the nonparametric test in statistics, this paper concluded that the best distribution of the mass of Saturn's regular moons is the lognormal distribution. A p-value close to one indicates that the best-fitting works well. According to the results of statistical inference, a mass function relationship that is different from the empirical Jeans' law was obtained, and on this basis, the dynamic equations of the variable mass restricted threebody problem were established. With the help of the Lindstedt-Poincaré perturbation method, the analytical expression of the approximate periodic orbit near the Lagrangian point L 3 was obtained. In addition, the effects of changes in scale parameter σ, the threebody interaction parameter k, and the variable-mass parameter γ on the periodic orbit and the position of the Lagrangian point L 3 were also discussed.
When σ increases, the Lagrangian point L 3 will gradually move away from the system's center-of-mass, and the degree of coincidence between the second-and third-order periodic orbits will also be significantly improved. In the range of its confidence interval [6.05109, 11.0738], the shape and amplitude of the orbit and the position of the Lagrangian point L 3 will remain almost unchanged, which shows that the statistical error has a small effect on the periodic orbit near the Lagrangian point L 3 . For the change of k, the difference between the second-and third-order periodic orbits will not change significantly, but the increase of the value of k will not only make the Lagrangian point L 3 move back to the center-of-mass of the system, but also, the amplitude of the orbit will increase. The closer γ is to one, the better the coincidence between the second-and third-order periodic orbits. In the mission of orbit design, it is believed that the second-order periodic orbit can replace the third-order periodic orbit as the initial approximation of the orbit design. Within a specific range, the change in γ will only affect the evolution of the periodic orbit and the Lagrangian point L 3 .