Abstract
Terms in the analytic expansion of the doubly averaged disturbing function for the circular restricted three-body problem using the Legendre polynomial are explicitly calculated up to the fourteenth order of semimajor axis ratio between perturbed and perturbing bodies in the inner case , and up to the fifteenth order in the outer case . The expansion outcome is compared with results from numerical quadrature on an equipotential surface. Comparison with direct numerical integration of equations of motion is also presented. Overall, the high-order analytic expansion of the doubly averaged disturbing function yields a result that agrees well with the numerical quadrature and with the numerical integration. Local extremums of the doubly averaged disturbing function are quantitatively reproduced by the high-order analytic expansion even when is large. Although the analytic expansion is not applicable in some circumstances such as when orbits of perturbed and perturbing bodies cross or when strong mean motion resonance is at work, our expansion result will be useful for analytically understanding the long-term dynamical behavior of perturbed bodies in circular restricted three-body systems.
1. Introduction
In the long tradition of celestial mechanics, the restricted three-body problem has occupied a fundamental role. In this problem, the mass of one of the three bodies is assumed to be small enough so that it does not affect the motion of the other two bodies. The restricted three-body problem is often considered on a rotating coordinate where central body and perturbing body are always located on the -axis. See Szebehely [1] for more general characteristics of the restricted three-body problem.
Among many variants of the restricted three-body problem, its circular version called the circular restricted three-body problem (hereafter referred to as CR3BP) has been studied particularly well, and it makes a basis for understanding solar system dynamics and many other fields in celestial mechanics. In this system, the perturbing body lies on a circular orbit around central body. As is well known, the degree of freedom of CR3BP becomes unity and the system turns into integrable once we average the disturbing function of the system by mean anomalies of perturbed and perturbing bodies. Theories of the so-called classical Lidov–Kozai cycle have been developed based on the integrable characteristics of the doubly averaged CR3BP [2–4], where stationary points of argument of pericenter around appear when the vertical component of the angular momentum of perturbed body is smaller than a certain value. The present paper deals with the doubly averaged CR3BP.
In the classical theory of the Lidov–Kozai cycle, the doubly averaged disturbing function is expanded using the Legendre polynomials of even orders. Putting , where is the semimajor axis of perturbed body and is that of perturbing body, only the lowest-order terms up to are considered in many of the studies along this line (e.g., [5–9]); it is the quadruple-order approximation. Recent studies of the so-called eccentric Lidov–Kozai mechanism (e.g., [10–13]) that deal with the eccentric restricted three-body problem (ER3BP), where the orbit of perturbing body has a finite eccentricity, are based on the octupole-order approximation of disturbing function up to . It is now getting better known that the inclusion of octupole terms in the disturbing function substantially changes the dynamical behavior of ER3BP. Even in CR3BP, the quadruple-order approximation is not accurate enough when is large. To eliminate this shortcoming, in the early days, Kozai [3] expanded the doubly averaged disturbing function of the inner CR3BP up to . More recently, Laskar and Boué [14] calculated analytic expansions of the general three-body disturbing function together with a practical method to compute the Hansen coefficients. Laskar and Boué [14] explicitly showed expressions of secular disturbing function up to for planar problems and up to for spatial problems.
In the present paper we will show specific expressions of the analytic expansion of the doubly averaged spatial disturbing function up to for the inner CR3BP and up to for the outer CR3BP using the Legendre polynomials. As most readers are aware, very wide varieties of studies have been already done on the analytic expansion of the doubly averaged disturbing function of CR3BP. Compared with previous literature, the present paper intentionally aims to be rather expository. The major purpose of this paper is to explicitly show expressions of the high-order analytic expansion of the doubly averaged disturbing function for CR3BP so that readers with interest can consult the high-order terms without going through algebraic manipulation by themselves. We also aim at analytically reproducing local extremums that secular disturbing function of CR3BP intrinsically has, particularly when or is large. Thus, even on a very basic subject like this, we have felt it advisable to give more details than would otherwise be necessary.
In Section 2 we give a brief description of the disturbing function of the three-body problem that we consider in the present paper. Section 3 goes to double averaging and analytic expansion of the disturbing function for CR3BP: general procedure (Section 3.1) and specific forms (from Sections 3.2 to 3.8). In Section 4 we show a comparison between the results obtained by the analytic expansion and by numerical quadrature. We also carried out direct numerical integration of equations of motion for comparison and show its result in Section 5. Section 6 is devoted to summary and discussion.
For readers’ convenience, before getting into the main sections let us quickly write down the basic equations of motion of the system that we deal with in the present paper. The differential equation that we will consider is the simple classical Newtonian equation of motionwhere indicates the position vector of the perturbed body and relates to the central mass. The disturbing function that plays a central role in this paper is denoted as . As for a literal definition of the doubly averaged disturbing function, particularly its direct part, we use the following one:where is related to the mass of perturbing body, and denote mean anomaly of perturbed and perturbing bodies, respectively, and is the osculating distance of the two bodies in space. Consult later sections for detailed definitions of the variables in the above equations. Needless to say, the considered system contains only three bodies.
2. Disturbing Function
In the present paper we categorize CR3BP in two cases: (i) the inner case where the orbit of the perturbed body is located inside that of a perturbing body , such as the Sun-asteroid-Jupiter system, and (ii) the outer case where the orbit of the perturbed body is located outside that of a perturbing body such as the Sun-Neptune-TNO system (TNO = Trans-Neptunian Object). The coorbital case is out of the scope of the present paper.
Following the long-term convention of celestial mechanics, in the present paper we express the disturbing function of CR3BP in relative coordinates where the origin is located on the primary body (Figure 1(a)). In this coordinate system, the disturbing function that describes the perturbation on the motion of an object with mass due to the motion of another mass has the following general form (e.g., [15, p. 228]):where with the gravitational constant , is the position vector of the mass with respect to the central body, is the position vector of the mass with respect to the central body, and . In what follows we will consider only the first term of the right-hand side of (3) which is often referred to as the direct part. The second term is called the indirect part. As is well known, the indirect part makes no contribution to long-term dynamics of the system because it vanishes after the double averaging procedure, unless nonnegligible mean motion resonances are at work and we cannot simply employ the double averaging procedure.
(a)
(b)
When designating as the angle between the vectors and , it is also well known that on the right-hand side of (3) can be expanded using the Legendre polynomials aswhen (i.e., the inner case), andwhen (i.e., the outer case). Once again, it is well known that the terms of and in (4) and (5) do not contribute to secular motion of the bodies, as they will vanish or become constant after the double averaging procedure. Hence in the remaining part of this paper we will just consider terms with in (4) and (5).
Readers find the expressions of the disturbing functions (3), (4), and (5) and their derivations in many textbooks such as Brouwer and Clemence [16], Danby [17], or Murray and Dermott [15]. In the inner case, the direct part of the disturbing function can be derived also in a more general way. Consider a general three-body system with three masses: primary , secondary , and tertiary . Now we use the Jacobi coordinate (e.g., [18]): measuring ’s position vector from , measuring ’s position vector from the barycenter of and , and is the angle between the vectors and . Naturally and have different origins, and the angle is different from in general (see Figure 1(b)). We assume . In this coordinate system, the equations of motion of and become (see [16, 19], for detailed derivation)whereare reduced masses used in the Jacobi coordinate system (e.g., [20, 21]), and is the common force functionwhereis the mass factor. Using the force function , this system can be written in a canonical form governed by a Hamiltonian. By expressing and as the semimajor axis of the orbits of the secondary and the tertiary (e.g., [22–25]), the Hamiltonian becomes
The first and the second terms of in (10) drive the Keplerian motion of the secondary and tertiary mass, respectively. Note that the third term which represents the mutual interaction of the secondary and the tertiary does not include terms of or . This is a typical consequence of the use of the Jacobi coordinate that separates the motions of three bodies into two separate binaries and their interactions by a single infinite series.
Now, let us think about the limit where is infinitesimally small; this would correspond to the inner R3BP. In this case, we divide the force function in (8) by in (7) for normalization by mass before taking the limit. Then the third term of becomes
Now we take the limit of , and the position of in Figure 1(b) gets overlapped with the position of . Thus we can replace for , for , for , and for and will end up with an expression equivalent to the direct part of the disturbing function of the inner case written in the relative coordinate (4).
On the other hand, deriving the form of the disturbing function of the outer case written in the relative coordinate (5) by taking the mass-less limit of Hamiltonian (10) is difficult, if not impossible. In the outer case, the origin of the position vector of the tertiary is the barycenter of the primary and secondary ( in Figure 1(b)). But would not be overlapped with the position of the primary regardless of the value of the tertiary’s mass , unless . Therefore, in the present paper, we will not mention the conversion of the general three-body Hamiltonian into the outer disturbing function written in the relative coordinate. In modern celestial mechanics, more and more methods for expanding disturbing function without using the relative coordinate are becoming available (e.g., [14, 26, 27]).
3. Doubly Averaged Disturbing Function for CR3BP
3.1. General Form
From (3) and (4), the direct part of the disturbing function for the inner CR3BP where becomes as follows:where we ignore the term of . In (12) we also ignore all the terms including the odd Legendre polynomials because they all vanish after the averaging procedure using mean anomaly of the perturbing body. Note that in the remaining part of this paper we will not consider the indirect part of the disturbing function either, as they do not have any secular dynamical contributions in nonresonant systems. Therefore we just use the variable for representing the entire part of the disturbing function.
Assuming there is no major resonant relationship between the mean motions of perturbed and perturbing bodies, we now try to get the double average of (12) over mean anomalies of both the bodies. Nonexistence of a resonant relationship means that the mean anomalies of perturbed and perturbing bodies (referred to as and in what follows) are independent from each other. The procedure to carry out double averaging of is straightforward as follows: Let us pick up the th term of and name it as . We have
First we average by mean anomaly of the perturbing body , aswhere
The angle is expressed by orbital angles through a relationship (e.g., [3, Eq. () in p. 592]) where , are true anomalies of the perturbed and perturbing bodies, , are arguments of pericenter of the perturbed and perturbing bodies, and is their mutual inclination measured at the node of the two orbits. We choose the orbital plane of the perturbing body as a reference plane for the entire system, and then and can be measured from the mutual node. Note that is not actually defined in CR3BP. Therefore, in (16) we regard as a single, fast-moving variable when we carry out averaging of (15). Practically, we can simply replace for in the discussion here.
To obtain of (15), we calculate the average of by asThen we average of (14) by mean anomaly of the perturbed body , as
If we switch the integration variable from to eccentric anomaly , (18) becomes
We can obtain the doubly averaged disturbing function for the outer CR3BP in the same way as above. From (3) and (5), the direct part of the disturbing function for the outer CR3BP becomes as follows:
Note that our definition of for the outer case (5), hence also in (20), may be different from what is seen in conventional textbooks (e.g., [15, Eq. () in p. 229]): Roles of the dashed quantities may be the opposite. This difference comes from the fact that conventional textbooks always assume , while we assume for the outer problem. This is because we make it a rule to always use dash for the quantities of perturbing body, whether it is located inside or outside the perturbed body.
Similar to the procedures that we went through for the inner CR3BP, we again assume that there is no major resonant relationship between mean motions of perturbed and perturbing bodies in the outer CR3BP. We then try to get the double average of over mean anomalies of both the bodies. Let us pick the th term of in (20) and name it . We have
First we average by mean anomaly of the perturbing body , aswhere is already defined in (15).
Then we average in (22) by mean anomaly of the perturbed body , as
If we switch the integration variable from to true anomaly , (23) becomes
Note that has the order of as in (23) and (24), not .
In Sections 3.2 to 3.8 we show the expressions of in (17), in (15), the Legendre polynomial in its original form, in (19), and in (24) for . We used Maple™ for algebraic manipulation to obtain the series of expressions. Note that in what follows we use instead of because it generally makes the formulas simpler. For this reason, some of the expressions look apparently different from what was presented in the previous literature in spite of their equivalence.
3.2.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
We just describe the resulting expressions of the expansion as follows: Let us emphasize again that the dashed quantities such as and are those of the perturbing body, whether its orbit is located inside or outside of the orbit of the perturbed body.
The expression in (28) shows the leading term of the doubly averaged disturbing function that causes the classical (circular) Lidov–Kozai cycle in the inner CR3BP that we have often seen in the previous literature. Meanwhile, the expression in (29) shows the leading term of the doubly averaged disturbing function for the outer CR3BP, but somehow we do not see it often. We should note that the leading term of the doubly averaged disturbing function for the outer problem (29) does not contain dependence on . The -dependence in the doubly averaged outer CR3BP first shows up in the next order: . This is why the Lidov–Kozai mechanism for the outer CR3BP is more subtle than the inner one, particularly when is small, and perhaps this is why we rarely see the expression in the literature.
3.3.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
Note that we now see the -dependence of the doubly averaged disturbing function for the outer CR3BP in the expression of (34).
3.4.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
3.5.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
3.6.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
3.7.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
3.8.
At this order, the corresponding component of the disturbing function for the inner problem is of , and that for the outer problem is of .
4. Comparison with Numerical Quadrature
To graphically show the validity of the high-order analytic expansion of the doubly averaged disturbing function that we have presented, let us carry out numerical quadrature for comparison. The literal definition of the doubly averaged disturbing function, particularly its direct part, becomes from (3)
Since the orbit of the perturbing body is circular in CR3BP, we can turn the double integral of (52) into a single integral using the complete elliptic integral of the first kind (e.g., [28–30]). To compare the numerical quadrature result with that from the analytic expansion of disturbing function, we draw equipotential curves on the plane. Since is an angle that intrinsically rotates, use of this polar-type coordinate is reasonable; therefore, it has often been used in the previous literature (e.g., [31–33]). All the numerical quadratures including the calculation of elliptic integral presented in this section are achieved using the functions implemented in GNU Scientific Library (GSL; http://www.gnu.org/software/gsl/).
As mentioned, the doubly averaged CR3BP makes an integrable system with just one degree of freedom. Then the state of the system is determined by just two constant parameters: (or ) and (square of) normalized vertical component of the perturbed body’s angular momentum defined as (e.g., [2, 34–37])
We selected several combinations of and drew equipotential curves using the analytically expanded doubly averaged disturbing function, as well as using the numerical quadrature (Figure 2). Note that from the definition of (53), the theoretical largest value of is achieved when as . Similarly, the theoretical smallest value of is achieved when as .
(a)
(b)
(c)
(d)
(e)
(f)
Figures 2(a) and 2(b) are for the inner case, and Figures 2(c)–2(f) are for the outer case. In each panel of Figure 2, the partial circle with red represents a line where the orbits of perturbed and perturbing bodies cross at the ascending node of the perturbed body. The partial circle with blue represents a line where the orbits of perturbed and perturbing bodies cross at the descending node of the perturbed body. These two circles are functions just of and (i.e., functions of ), and the relationship between , , , and is expressed as (e.g., [37–40])In the denominator of the right-hand term of (54), the positive sign takes place when the orbit intersect happens at the ascending node of perturbed body (when , where denotes perturbed body’s true anomaly). The negative sign takes place when the encounter happens at the descending node of perturbed body (when ). Using the variables and , (54) can be rewritten asand it is more obvious that (55) represents a pair of circles on the plane.
Note that the orbit-crossing lines introduce singularities in -body Hamiltonian (e.g., [41, 42]), and (or ) is continuous but not regular across the orbit-crossing lines. Therefore the analytic expansion of the doubly averaged disturbing function is not applicable to the region very close to the orbit-crossing lines (e.g., [43]). Note also that each panel in Figure 2 has an individual contour interval, which represents a fixed interval of disturbing potential in each of the systems. Therefore, regions with dense contours imply that the potential gradient is steep there.
Since the high-order analytic expansion of disturbing function must have its significance when (or ) is large, we chose relatively large (or ) so that we cannot ignore their high-order powers. For example, when and , which we may ignore. But when and , , which is not ignorable. In Figure 2(a) where , effect of the high-order analytic expansion is well exhibited. As is widely known, there is no stationary point for argument of pericenter while in the doubly averaged inner CR3BP at the quadruple-order approximation . However when becomes large, acquires stationary points at even when . This is what we see in the quadrature result in the leftmost panel in Figure 2(a). Kozai [3] derived expressions of high-order analytic expansion of the doubly averaged disturbing function for the inner CR3BP up to . But the stationary points of seen in the rightmost panel of Figure 2(a) (for quadrature) are not reproduced even by using the analytic terms. We confirmed that the stationary points at in this parameter set first appear in the analytic approximation up to , and the value of at the stationary points is better reproduced using the analytic expansion of , as seen in Figure 2(a). Note that other sets of stationary points at seen in the quadrature result in the leftmost panel are not reproduced by the analytic expansion of disturbing function, since they are out of the orbit-crossing lines.
Note that we did not draw all the contours near the outer boundary indicated by the dashed circles, particularly in the panels in the left three columns. This is mainly because the data size of the figure would be too large if we drew all the contours with a fixed interval of disturbing potential toward the outer boundary, and also because the contour intervals would become too narrow. Moreover, we believe that we can grab major topological patterns of disturbing potential in each panel even if we draw contours only partially, at least inside the orbit-crossing lines.
Figure 2(b) is for . They exemplify some limitations of the high-order analytic expansion of disturbing function that we have carried out. Although a saddle point seen in the quadrature result at the origin is reproduced in the results of the analytic expansion of disturbing function, we see spurious local extremums along the -axis near the boundary of in the case. The spurious local extremums are not so remarkable in the case (although they surely exist), but there seem to exist spurious saddle points around that do not appear in the quadrature result. Thus we say that the analytic expansion of the doubly averaged disturbing function for the inner CR3BP, at least up to these orders, is applicable to systems with moderate eccentricity when is large. Note that it apparently seems that the leftmost panel for the analytic expansion of the order of looks the closest to the quadrature result. However, we think this is just a coincidence, because the expansion includes only the terms of as in (28), and would not be influenced by terms of with .
Figure 2(c) is for the doubly averaged outer CR3BP when . In this case, the high-order analytic expansion of the doubly averaged disturbing function reproduces the quadrature result inside the orbit-crossing lines. However, the medium-order expansion yields artificial saddle points along the -axis. We confirmed that these spurious saddle points show up in the expansion of too (panels are not shown here). We thus say that the highest-order analytic expansion of pushes these saddle points out of the orbit-crossing lines, and the equipotential curves become similar to the quadrature result, at least inside the orbit-crossing lines where is moderate.
In Figure 2(d) when , we see three local extremums in the quadrature result: a pair of local minima along the -axis together with a saddle point at the origin . These local extremums cannot be reproduced by the low-order analytic expansion such as or ; the local extremums show up only in the expansion of or higher, although the locations of the local minima along the -axis are slightly different from the quadrature result. When becomes even smaller, the difference in the topology of the equipotential curves increases, as seen in Figure 2(e) where . A pair of local minima along the -axis persists, but the origin becomes a local maximum and a pair of saddle points takes place along the -axis. We confirmed that this topology is only reproduced when using the analytic expansion of disturbing function at our highest-order , although the locations of the local extremums are slightly different from the quadrature result.
When gets even larger and is in a certain range of moderate values, a pair of local minima show up along the -axis inside the orbit-crossing lines together with a saddle point at the origin . This is the case of Figure 2(f) when . The topology of the local extremums is reproduced by the analytic expansion of disturbing function at our highest-order . However, we should note a point here. In the quadrature result of the rightmost panel in Figure 2(f), the pair of local minima along the -axis appear when the value of is slightly different from the analytic result . We have not yet figured out the cause of the discrepancy of the values that are required for the local minima to appear in the analytic expansion result and in the numerical quadrature result. It is possible that even higher-order analytic expansions of disturbing function may diminish the discrepancy.
As a summary of this section, we say that the high-order analytic expansion of the doubly averaged disturbing function of CR3BP that has been presented in the present paper overall seems valid inside the orbit-crossing lines even when (or ) is large, as long as eccentricity of perturbed body is moderate. Under these conditions, local extremums of the doubly averaged disturbing function for CR3BP are quantitatively reproduced by the high-order analytic expansion.
5. Comparison with Numerical Integration
In Section 4 we inspected the validity of the high-order analytic expansion of the doubly averaged disturbing function by comparing its outcome with that from numerical quadrature. The numerical quadrature defined by the double integral (52) yields a numerical proxy of an integrable CR3BP system with one degree of freedom, and it is safe for us to directly compare the quadrature result with that from the analytic expansion of disturbing function. However, equations of motion of three-body systems are generally not integrable, having short-term oscillations and various mean motion resonances intact without being averaged out. To demonstrate the practical usefulness of our analytic expansion of the doubly averaged disturbing function, we carried out direct numerical integration of equations of motion using the same parameter combinations (or ) and as in the previous section. Based on the comparison, we will again give considerations on the validity of the high-order analytic expansion of the doubly averaged disturbing function for CR3BP.
5.1. Model and Method
We consider motion of small solar system bodies with infinitesimally small mass (asteroid, comet, etc.) undergoing gravitational perturbation from a major planet on a circular orbit. Using the same six sets of parameters in the framework of CR3BP as in Figure 2, we calculate orbit propagation of the small bodies by numerical integration for ten to two hundred million years. The parameter sets are tabulated in Table 1. The main result of the numerical integration is presented and discussed in Section 5.2.
In Table 1, mass ratios between perturbing body and central body for (a), (b), and (c) are close to the ratio of the Earth-Moon total mass and the solar mass (~ of the mass ratio between Jupiter and the Sun). The value for (f) is nearly of the mass ratio of Jupiter and the Sun. The reason for these small values is to avoid occurrence of orbital instability of the small body due to close encounters with the perturbing planet near the orbit-crossing lines. This is necessary because and are large in these four cases. We also chose small values for the stepsize of numerical integration for these cases by the same reason. As for (d) and (e), the value is close to of the mass ratio of Jupiter and the Sun. Since is relatively small in these two cases, the probability of close encounters between perturbed and perturbing bodies is lower, and we can choose a larger mass ratio as well as larger stepsize . Here, readers should recall the fact that the influence of the mass ratio is limited to the timescale of orbital evolution of perturbed body, and it does not affect the topology of relative orbit, particularly in the doubly averaged CR3BP. It is obvious from the function form of the doubly averaged disturbing function (19) or (24) where the perturber’s mass serves just as a constant coefficient through . Therefore, the trajectory shape of the perturbed body remains the same even if we change the mass of the perturber from to (where is a positive real number); only the timescale of orbital evolution would change from to . This is why we can arbitrarily change the mass ratio in numerical integration as long as what we need to know is the secular topology of trajectories of the perturbed body, not its time evolution sequence.
In the numerical integration, the semimajor axis of the perturbing body AU is common to all the cases, which is a substitute for Jupiter’s semimajor axis. The perturber’s mean anomaly (practically equivalent to mean longitude) is set to zero at time . By the definition of CR3BP, the perturber’s eccentricity and inclination are both set to zero. As for the perturbed body’s initial orbital elements, semimajor axis is automatically determined from the value of or . We set the initial values of eccentricity and argument of pericenter of perturbed body along the - and -axis on the plane: we selected of perturbed body from to with the interval of 0.01. As for , we used four values of , , , and . The values of the perturbed body’s are automatically determined by relation (53). Initial values of the perturbed body’s longitude of ascending node are set to 0, and those for its initial mean anomaly are set to .
The differential equation to be solved is the simple classical Newtonian equation of motion of the perturbed body. For the inner case, it iswith . For the outer case we just replace for . As for the numerical integration scheme, we use the regularized mixed-variable symplectic method [44] based on the so-called Wisdom–Holman map [20, 45]. In this scheme the stepsize of numerical integration is automatically reduced when a close encounter between two bodies happens, trying to keep the calculation accuracy high. The values listed in Table 1 are the nominal stepsizes used when no close encounter takes place between bodies. We have implemented this scheme based on the code SWIFT [44] and used it in our previous studies where close encounters between asteroids and terrestrial planets take place very often [46, 47], as well as where almost no close encounter happens [48, 49]. We have confirmed that our implementation of this scheme does not have a practical problem in accuracy or efficiency when applied to these dynamical systems. We adopt the Gaussian units where the Gaussian gravitational constant , the unit of time is an ephemeris day, and the unit of length is the astronomical unit (AU), so that the mass of the primary body becomes unity (e.g., [16, 50]).
5.2. Integration Result
We show the main results of our numerical integration in Figure 3 whose parameters are described in Table 1. The perturbed body’s values are plotted on the panels by black and green dots (see below for the reason for this distinction) with the fixed output interval designated in Table 1. We chose empirically so that the output data amount does not get too large, while avoiding any kind of aliasing in the plots. The following settings are all common to both Figures 2 and 3: the coordinate system , the ranges of - and -axis, meaning of the black dashed circles, and that of the red and blue partial circles.
When inspecting Figure 3, we would like readers to notice two points: One is that we drew trajectories of perturbed bodies in black only when their orbital motion remains stable from the beginning to the end of the integration period . This is to facilitate the comparison between the quadrature result and the integration result. Here we define the “stable” orbital motion if the semimajor axis does not exhibit any major changes throughout the integration period . Otherwise we categorize the orbital motion as “unstable.” This categorization is done by the naked eye, plotting the time sequence of for all the perturbed bodies used in the integration (320 bodies per panel, 1920 bodies in total). Browsing through all the plots of , we subjectively distinguish “stable” and “unstable” orbits. Hence there can be some ambiguity in the categorization. However, let us say that the difference between stable and unstable orbits is rather clear in Figure 3, and the subjective distinction between the black and the green trajectories does not affect the discussion here.
The second point is that we did not draw some of the trajectories that are largely overlapped with others. This is to prevent the panels from becoming too busy with too many similar trajectories, as well as to keep the figure file size reasonably small. As we mentioned before, starting points of the numerical integration are aligned on the - and -axis as with a fixed interval of eccentricity. Due to the intrinsic symmetry that the disturbing function has, different initial condition can generate almost identical trajectories. For example, in Figure 3(c) the initial conditions and generate almost identical trajectories. For eliminating redundancy in figure and reducing the file size, we did not draw the trajectory from the initial condition and just drew that from . We carried out this kind of selection of overlapped trajectories with great care so that our subjective selection never alters general topological patterns in the panels.
Comparing Figures 3 and 2, we say that their agreement is overall good. In Figure 3(a) for the inner CR3BP, existence of stationary points at is confirmed by the numerical integration. The eccentricity value of the stationary points is close to what the analytic expansion and the numerical quadrature had yielded in Figure 2(a). Existence of stationary points at outside the orbit-crossing lines is also confirmed, although most of the trajectories near the orbit-crossing lines become unstable due to close encounters with the perturbing body.
Similarly in Figure 3(b), existence of stationary points at as well as those at is confirmed by the numerical integration. There must be a saddle point at the origin , but most of the trajectories near the origin are unstable in Figure 3(b) due to the proximity to the orbit-crossing lines, and it is not easy to visually confirm it.
In Figure 3(c) for the outer CR3BP, existence of a stationary point at the origin as well as those at is confirmed by the numerical integration. It is clearly seen that trajectories can easily become unstable around the two intersections of the red and blue orbit-crossing lines along the -axis.
In Figures 3(d) and 3(e), the existence of stationary points is again confirmed by the numerical integration. We see a pair of stationary points at in both panels. The origin in Figure 3(d) makes a saddle point, as suggested by the quadrature (Figure 2(d)). The origin in Figure 3(e) makes a local maximum, which is also consistent with the quadrature result (Figure 2(e)). A pair of saddle points seen along the -axis in Figure 3(e) has also been suggested by the quadrature (Figure 2(e)). Although the location of the saddle points seen in Figure 3(e) is slightly closer to the origin compared with that in Figure 2(e), we would say they are consistent. In both the panels of Figures 3(d) and 3(e), trajectories get unstable and the green dots are scattered as they approach the orbit-crossing lines.
As we have seen in Figures 3(a)–3(e), the result from the direct numerical integration of the equations of motion overall agrees well with the result from the numerical quadrature, as well as with that from the analytic expansion of the doubly averaged disturbing function. We thus say that this comparison largely justifies the validity of use of the analytic expansion of the doubly averaged disturbing function for CR3BP presented in Section 3.
However, in Figure 3(f) we encounter with a typical system where the use of doubly averaged disturbing function, either through analytic expansion or quadrature, is not appropriate. As is clearly seen, the topological pattern of trajectories in Figure 3(f) is very different from the equipotential contours shown in Figure 2(f). A pair of stationary points shows up at in Figure 3(f) that we do not see in Figure 2(f). Also, although we see stationary points along the -axis in Figure 3(f), the shape of the trajectories is fairly different from what is observed in Figure 2(f). We then plotted the whole phase space for this system in Figure 4. As we see, the trajectories outside the orbit-crossing lines are very similar between Figure 4(a) (quadrature) and Figure 4(b) (integration). But the complicated trajectory pattern inside the orbit-crossing lines takes place only in the numerical integration result (Figure 4(b)).
(a)
(b)
So far we have not completely figured out the cause of the difference in trajectory patterns between Figures 2(f) and 3(f). We have confirmed that there is no problem in accuracy of the numerical integration that is produced in Figure 3. As far as we have been able to investigate, the apparently peculiar trajectories seen in Figure 3(f) are possibly related to the existence of a mean motion resonance. We selected some of the periodic orbits observed inside the orbit-crossing lines in Figure 4(b) and plotted them in separate panels in Figure 5 (the panels in its left column). For each of the periodic orbits, we calculated an argument , where is mean longitude of the perturbing body and is that of the perturbed body. As a result, we found that in the top two trajectories (Figures 5(a) and 5(b)) the argument librates around , although the libration amplitude is relatively large (see the panels in the middle column of Figure 5). This indicates that the systems in Figures 5(a) and 5(b) are trapped in the 5 : 7 mean motion resonance. System (c) seems close to this resonance too. It is well known that the existence of strong mean motion resonance can substantially change the shape of a disturbing potential, together with the location of local extremums, compared with nonresonant cases (e.g., [51–53]). To deal with resonant systems like Figure 3(f) rigorously, we would need to employ different methods that are particular to each of the resonant systems (e.g., [54, 55]).
(a)
(b)
(c)
(d)
(e)
Although all the five perturbed bodies whose trajectories are shown in Figure 5 started with the same initial semimajor axis (), we see a systematic difference in the average value of the semimajor axis between 5(a–c) and 5(d, e) (see the panels in the right column of Figure 5). The short-term oscillation of with relatively large amplitude seen in Figures 5(a), 5(b), and 5(c) is probably related to the 5 : 7 mean motion resonance. It is possible that the oscillation of seen in Figures 5(d) and 5(e) is also related to this resonance or others that we have not yet found.
We will pursue the above problem in our forthcoming publications. Several possibilities might account for the peculiar trajectories seen in Figure 4(b), such as the effect of high-order mean motion resonance (e.g., [56]) and its combination with the Lidov–Kozai libration (e.g., [57]). A kind of secondary resonance may be embedded in this system (e.g., [58]). We should be aware that more and more mean motion resonances show up as or approaches 1 (e.g., [59]), and treatment of such systems would need greater care.
For readers’ interest, we also plotted trajectories of the perturbed body where is its inclination and is its longitude of ascending node. The six panels from (a) to (f) in Figure 6 correspond to each in Figure 3 with the same color pattern (stable orbits are drawn in black, and unstable orbits are drawn in green). Recalling the fact that the theoretical smallest value of takes place when as (see Section 4), it is clear that the theoretical largest value of is when (note that this value happens to be the same as ). We drew the outer boundary by the black dashed circles in Figure 6.
In the doubly averaged CR3BP system that we now consider, eccentricity and inclination are correlated with each other through the definition of in (53). However, since longitude of ascending node is not included in the disturbing function does not have a particular correlation with eccentricity , inclination , or argument of pericenter . The panels of Figures 6(b) and 6(c) typically exemplify this circumstance where we see the trajectories of perturbed bodies which are not closed on the plane; and are not correlated, and these trajectories do not compose equipotential contours.
For more demonstrations of the time variation of and , we picked one body per each of the panels from Figure 6(a) to Figure 6(f) and showed the time variation of their , , , and in Figure 7. Initial values of of each of the bodies are written in the caption of Figure 7. Our choice of the particular body among the 320 bodies in each of the panels is almost arbitrary: we basically chose the bodies whose argument of pericenter librates around or 0 (a, b, c, d, and f), but we chose a body whose argument of pericenter circulates from 0 to (e).
(a)
(b)
(c)
(d)
(e)
(f)
The panels in Figure 7 immediately tell us the correlation between eccentricity and inclination , which is an obvious outcome of the preservation of in (53). The argument of pericenter is correlated to them, as we have seen in the study of the doubly averaged disturbing function in Section 3. The rightmost panels show the fact that only longitude of ascending node has no correlation with other elements; it shows a monotonic decrease in all the panels.
We do not show the analytic counterpart of Figure 6 or Figure 7 because they are out of the scope of the present paper. Readers with particular interest may consult classic literature such as Moiseev [34, 60] or Lidov [2] for how to formally obtain analytic solutions of time evolution of . For example, in Lidov [2], the author uses the mutual dependence between , , and and constructs a definite integral for quadrature that formally depicts the function form of depending on time. Incidentally, if we limit ourselves to the doubly averaged inner CR3BP at the quadruple-order analytic approximation (; see Section 3.2), analytic solution of and that of (hence those of and ) using elliptic functions have already been derived in much of the literature (e.g., [3, 5–7, 61, 62]).
6. Summary and Discussion
Through the comparison of the results obtained from analytic expansion of the doubly averaged disturbing function, numerical quadrature, and direct numerical integration of the equations of motion, overall we validated the use of our high-order analytic expansion of the doubly averaged disturbing function both for the inner and outer CR3BP, even when the semimajor axis ratio (or ) is large. Our result tells us that the high-order analytic expansion is not applicable in several circumstances: For example, when a perturbed body’s eccentricity approaches its theoretical largest value, when orbits of perturbed and perturbing bodies cross, and when strong mean motion resonance is at work. Particularly, occurrence of strong mean motion resonance can significantly change global topology of disturbing potential, and the doubly averaged method itself goes out of use. Nevertheless, we think our original objective to analytically exemplify dynamical characteristics of doubly averaged disturbing function for CR3BP in this paper is largely fulfilled. We thus conclude that the high-order analytic expansion of the doubly averaged disturbing function is suitable for expository purposes to account for the dynamical characteristics of secular CR3BP even when (or ) is large, as long as the system is free from serious mean motion resonance.
Let us just mention a point regarding one limitation of the high-order analytic expansion of the doubly averaged disturbing function. As seen in Sections 4 and 5, sometimes the locations of local extremums of the analytically expanded doubly averaged disturbing function differ from those obtained by the numerical quadrature and numerical integration, particularly in the outer problem such as in Figure 2(d) or 2(e). In the outer case, the factor in (24) substantially enhances the magnitude of high-order terms in the analytically expanded disturbing function when eccentricity is large. For example, when the factor becomes 8.65 , 27.3 , and 48.6 . When , the factor reaches 2127 , 126626 , and 977049 , possibly causing artificial oscillation of the higher-order analytic terms in the large- region.
Incidentally, we would like readers to be aware of a recent study that revealed that solutions of the three-body problem obtained from the double averaging procedure can substantially deviate from true solution due to accumulation of the effect of short-term oscillation [63]. This can happen particularly when the mass of the perturbing body is not negligibly small compared to the primary mass.
A Maple script that yields our analytic expansion in the present paper (see (25)–(51)) is available from the author upon request.
Competing Interests
The author declares no competing financial interests associated with this article.
Acknowledgments
The author has benefited from stimulating enlightenment through discussions with Katsuhito Ohtsuka, Arika Higuchi, Akihiko Fujii, Ruslan Salyamov, Fumi Yoshida, and Renu Malhotra. Detailed and constructive review by Yolande McLean has considerably improved the English presentation of this paper. Algebraic manipulations with Maple, numerical quadrature operations, and orbit propagation by numerical integration carried out in the present paper were all performed at Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan (NAOJ). This study is supported by a pair of JSPS Kakenhi Grants, JP25400458/2013–2016 and JP16K05546/2016–2018.