# High-Order Analytic Expansion of Disturbing Function for Doubly Averaged Circular Restricted Three-Body Problem

Authors: Takashi Ito
Journal: Advances in Astronomy
 (2016)

---

## 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 (α<1), and up to the fifteenth order in the outer case (α>1). 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.


---

## Body

## 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 thex-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 g around ±π/2 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α=a/a′, where a is the semimajor axis of perturbed body and a′ is that of perturbing body, only the lowest-order terms up to O(α2) 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 O(α3). 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 O(α8). 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 O(α14) for planar problems and up to O(α5) for spatial problems.In the present paper we will show specific expressions of the analytic expansion of the doubly averaged spatial disturbing function up toO(α14) for the inner CR3BP (α<1) and up to O(α′15) for the outer CR3BP (α′=α-1<1) 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 Section2 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 motion(1)d2rdt2+μrr3=∇R,where r 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 R. As for a literal definition of the doubly averaged disturbing function, particularly its direct part, we use the following one:(2)Rl′l=μ′4π2∬02πdldl′Δ,where μ′ is related to the mass of perturbing body, l and l′ 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(α<1), 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 (α′=α-1<1) such as the Sun-Neptune-TNO system (TNO = Trans-Neptunian Object). The coorbital case (α=1) 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 functionR 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 m due to the motion of another mass m′ has the following general form (e.g., [15, p. 228]):(3)R=μ′Δ-μ′r·r′r′3,where μ′=Gm′ with the gravitational constant G, r is the position vector of the mass m with respect to the central body, r′ is the position vector of the mass m′ with respect to the central body, and Δ=r′-r. 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.Figure 1
A schematic illustration of the three-body configuration considered here in two kinds of coordinate. (a) Relative coordinate. Both the vectorsr and r′ originates from the central mass denoted as P. (b) Jacobi coordinate. The vector r1 is originated from the primary mass (m0), and the vector r2 is originated from the barycenter of the primary mass and the secondary mass (m1) denoted as G.
(a) (b)When designatingS as the angle between the vectors r and r′, it is also well known that 1/Δ on the right-hand side of (3) can be expanded using the Legendre polynomials Pj as(4)1Δ=1r′∑j=0∞rr′jPjcos⁡S,when r<r′ (i.e., the inner case), and(5)1Δ=1r∑j=0∞r′rjPjcos⁡S,when r>r′ (i.e., the outer case). Once again, it is well known that the terms of j=0 and j=1 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 j≥2 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 m0, secondary m1, and tertiary m2. Now we use the Jacobi coordinate (e.g., [18]): measuring m1’s position vector r1 from m0, measuring m2’s position vector r2 from the barycenter of m0 and m1, and S12 is the angle between the vectors r1 and r2. Naturally r1 and r2 have different origins, and the angle S12 is different from S in general (see Figure 1(b)). We assume r1<r2. In this coordinate system, the equations of motion of m1 and m2 become (see [16, 19], for detailed derivation)(6)m~1d2r1dt2=∂F∂r1,m~2d2r2dt2=∂F∂r2,where(7)m~1=m0m1m0+m1,m~2=m0+m1m2m0+m1+m2are reduced masses used in the Jacobi coordinate system (e.g., [20, 21]), and F is the common force function(8)F=Gm0m1r1+m0+m1m2r2+1r2∑j=2∞Mjr1r2jPjcos⁡S12,where(9)Mj=m0m1m2m0j-1--m1j-1m0+m1jis the mass factor. Using the force function F, this system can be written in a canonical form governed by a Hamiltonian. By expressing a1 and a2 as the semimajor axis of the orbits of the secondary and the tertiary (e.g., [22–25]), the Hamiltonian H becomes(10)H=Gm0m12a1+Gm0+m1m22a2+Gr2∑j=2∞Mjr1r2jPjcos⁡S12.The first and the second terms ofH 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 j=0 or j=1. 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 wherem1 is infinitesimally small; this would correspond to the inner R3BP. In this case, we divide the force function F in (8) by m~1 in (7) for normalization by mass before taking the limit. Then the third term of F becomes(11)F3rdm~1=Gm2r2∑j=2∞m0j-1--m1j-1m0+m1j-1r1r2jPjcos⁡S12.Now we take the limit ofm1→0, and the position of G in Figure 1(b) gets overlapped with the position of m0. Thus we can replace r1 for r, r2 for r′, S12 for S, and m2 for m′ and will end up with an expression equivalent to the direct part of the disturbing function R 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 (m2→0) of Hamiltonian (10) is difficult, if not impossible. In the outer case, the origin of the position vector of the tertiary (r2) is the barycenter of the primary and secondary (G in Figure 1(b)). But G would not be overlapped with the position of the primary regardless of the value of the tertiary’s mass m2, unless m1→0. 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 r′=a′ becomes as follows:(12)R=μ′a′∑n=1∞ra′2nP2ncos⁡S,where we ignore the term of n=0. In (12) we also ignore all the terms including the odd Legendre polynomials (P1,P3,P5,P7,P9,…) 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 R 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 ofR (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 l and l′ in what follows) are independent from each other. The procedure to carry out double averaging of R is straightforward as follows: Let us pick up the nth term of R and name it as R2n. We have(13)R2n=μ′a′ra′2nP2ncos⁡S.First we averageR2n by mean anomaly of the perturbing body l′, as(14)R2nl′=μ′a′ra′2nP2nl′,where(15)P2nl′=12π∫02πP2ncos⁡Sdl′.The angleS is expressed by orbital angles through a relationship (e.g., [3, Eq. (7) in p. 592]) (16)cos⁡S=cos⁡f+gcos⁡f′+g′+cos⁡isin⁡f+gsin⁡f′+g′,where f, f′ are true anomalies of the perturbed and perturbing bodies, g, g′ are arguments of pericenter of the perturbed and perturbing bodies, and i 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 g and g′ can be measured from the mutual node. Note that g′ is not actually defined in CR3BP. Therefore, in (16) we regard f′+g′ as a single, fast-moving variable when we carry out averaging of (15). Practically, we can simply replace ∫dl′ for ∫df′ in the discussion here.To obtainP2nl′ of (15), we calculate the average of cos2n⁡S by l′ as(17)cos2n⁡Sl′=12π∫02πcos2n⁡Sdl′.Then we average R2nl′ of (14) by mean anomaly of the perturbed body l, as(18)R2nl′l=μ′a′aa′2n12π∫02πra2nP2nl′dl.If we switch the integration variable froml to eccentric anomaly u, (18) becomes(19)R2nl′l=μ′a′aa′2n12π∫02π1-ecos⁡u2n+1P2nl′du.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 R′ for the outer CR3BP becomes as follows:(20)R′=μ′r∑n=1∞a′r2nP2ncos⁡S.Note that our definition of1/Δ for the outer case (5), hence also in (20), may be different from what is seen in conventional textbooks (e.g., [15, Eq. (6.22) in p. 229]): Roles of the dashed quantities (r′,μ′) may be the opposite. This difference comes from the fact that conventional textbooks always assume r/r′<1, while we assume r/r′>1 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 ofR′ over mean anomalies of both the bodies. Let us pick the nth term of R′ in (20) and name it R2n′. We have(21)R2n′=μ′ra′r2nP2ncos⁡S.First we averageR2n′ by mean anomaly of the perturbing body l′, as(22)R2n′l′=μ′ra′r2nP2nl′,where P2nl′ is already defined in (15).Then we averageR2n′l′ in (22) by mean anomaly of the perturbed body l, as(23)R2n′l′l=μ′a′a′a2n+112π∫02πar2n+1P2nl′dl.If we switch the integration variable froml to true anomaly f, (23) becomes(24)R2n′l′l=μ′a′a′a2n+11-e2-2n+1/22π∫02π1+ecos⁡f2n-1P2nl′df.Note thatR2n′l′l has the order of O(α′2n+1) as in (23) and (24), not O(α′2n).In Sections3.2 to 3.8 we show the expressions of cos2n⁡Sl′ in (17), P2nl′ in (15), the Legendre polynomial P2n(x) in its original form, R2nl′l in (19), and R2n′l′l in (24) for 2n=2,4,6,8,10,12,14. We used Maple™ for algebraic manipulation to obtain the series of expressions. Note that in what follows we use sin⁡i instead of cos⁡i 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.2n=2

At this order, the corresponding component of the disturbing function for the inner problemR2 is of O(α2), and that for the outer problem R2′ is of O(α′3).We just describe the resulting expressions of the expansion as follows: Let us emphasize again that the dashed quantities such asμ′ and a′ are those of the perturbing body, whether its orbit is located inside or outside of the orbit of the perturbed body.(25)cos2⁡Sl′=14sin2⁡icos⁡2f+g-14sin2⁡i+12,(26)P2x=32x2-12,(27)P2l′=38sin2⁡icos⁡2f+g-38sin2⁡i+14,(28)R2l′l=μ′a′aa′21516e2sin2⁡icos⁡2g-1163e2+23sin2⁡i-2,(29)R2′l′l=μ′a′a′a31-e2-3/2-38sin2⁡i+14.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 g. The g-dependence in the doubly averaged outer CR3BP first shows up in the next order: 2n=4. 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.2n=4

At this order, the corresponding component of the disturbing function for the inner problemR4 is of O(α4), and that for the outer problem R4′ is of O(α′5).(30)cos4⁡Sl′=364sin4⁡icos⁡4f+g-316sin2⁡isin2⁡i-2cos⁡2f+g+964sin4⁡i-38sin2⁡i+38,(31)P4x=358x4-154x2+38,(32)P4l′=105512sin4⁡icos⁡4f+g-15128sin2⁡i7sin2⁡i-6cos⁡2f+g+315512sin4⁡i-4564sin2⁡i+964,(33)R4l′l=μ′a′aa′466154096e4sin4⁡icos⁡4g-3151024e2sin2⁡ie2+27sin2⁡i-6cos⁡2g+9409615e4+40e2+835sin4⁡i+40sin2⁡i+8,(34)R4′l′l=μ′a′a′a51-e2-7/2-45512e2sin2⁡i7sin2⁡i-6cos⁡2g+910243e2+235sin4⁡i-40sin2⁡i+8.Note that we now see theg-dependence of the doubly averaged disturbing function for the outer CR3BP in the expression of (34).


### 3.4.2n=6

At this order, the corresponding component of the disturbing function for the inner problemR6 is of O(α6), and that for the outer problem R6′ is of O(α′7).(35)cos6⁡Sl′=5512sin6⁡icos⁡6f+g-15256sin4⁡isin2⁡i-2cos⁡4f+g+15512sin2⁡i5sin4⁡i-16sin2⁡i+16cos⁡2f+g-5256sin2⁡i-25sin4⁡i-8sin2⁡i+8,(36)P6x=23116x6-31516x4+10516x2-516,(37)P6l′=11558192sin6⁡icos⁡6f+g-3154096sin4⁡i11sin2⁡i-10cos⁡4f+g+5258192sin2⁡i33sin4⁡i-48sin2⁡i+16cos⁡2f+g-57754096sin6⁡i+47252048sin4⁡i-525512sin2⁡i+25256,(38)R6l′l=μ′a′aa′6495495131072e6sin6⁡icos⁡6g-1039565536e4sin4⁡i3e2+1011sin2⁡i-10cos⁡4g+1575131072e2sin2⁡i15e4+80e2+4833sin4⁡i-48sin2⁡i+16cos⁡2g-256553635e6+210e4+168e2+16231sin6⁡i-378sin4⁡i+168sin2⁡i-16,(39)R6′l′l=μ′a′a′a71-e2-11/2-157565536e4sin4⁡i11sin2⁡i-10cos⁡4g+262532768e2sin2⁡ie2+233sin4⁡i-48sin2⁡i+16cos⁡2g-253276815e4+40e2+8231sin6⁡i-378sin4⁡i+168sin2⁡i-16.


### 3.5.2n=8

At this order, the corresponding component of the disturbing function for the inner problemR8 is of O(α8), and that for the outer problem R8′ is of O(α′9).(40)cos8⁡Sl′=3516384sin8⁡icos⁡8f+g-352048sin6⁡isin2⁡i-2cos⁡6f+g+354096sin4⁡i7sin4⁡i-24sin2⁡i+24cos⁡4f+g-352048sin2⁡isin2⁡i-27sin4⁡i-16sin2⁡i+16cos⁡2f+g+122516384sin8⁡i-175512sin6⁡i+315512sin4⁡i-3564sin2⁡i+35128,(41)P8x=6435128x8-300332x6+346564x4-31532x2+35128,(42)P8l′=2252252097152sin8⁡icos⁡8f+g-15015262144sin6⁡i15sin2⁡i-14cos⁡6f+g+24255524288sin4⁡i65sin4⁡i-104sin2⁡i+40cos⁡4f+g-11025262144sin2⁡i143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2f+g+78828752097152sin8⁡i-52552565536sin6⁡i+36382565536sin4⁡i-110258192sin2⁡i+122516384,(43)R8l′l=μ′a′aa′82737609875268435456e8sin8⁡icos⁡8g-1073572533554432e6sin6⁡i3e2+1415sin2⁡i-14cos⁡6g+1734232567108864e4sin4⁡ie4+8e2+865sin4⁡i-104sin2⁡i+40cos⁡4g-60637533554432e2sin2⁡i7e6+70e4+112e2+32143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2g+1225268435456315e8+3360e6+6048e4+2304e2+1286435sin8⁡i-13728sin6⁡i+9504sin4⁡i-2304sin2⁡i+128,(44)R8′l′l=μ′a′a′a91-e2-15/2-10510516777216e6sin6⁡i15sin2⁡i-14cos⁡6g+16978516777216e4sin4⁡i3e2+1065sin4⁡i-104sin2⁡i+40cos⁡4g-7717516777216e2sin2⁡i15e4+80e2+48143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2g+12253355443235e6+210e4+168e2+166435sin8⁡i-13728sin6⁡i+9504sin4⁡i-2304sin2⁡i+128.


### 3.6.2n=10

At this order, the corresponding component of the disturbing function for the inner problemR10 is of O(α10), and that for the outer problem R10′ is of O(α′11).(45)cos10⁡Sl′=63131072sin10⁡icos⁡10f+g-31565536sin8⁡isin2⁡i-2cos⁡8f+g+315131072sin6⁡i9sin4⁡i-32sin2⁡i+32cos⁡6f+g-31516384sin4⁡isin2⁡i-23sin4⁡i-8sin2⁡i+8cos⁡4f+g+31565536sin2⁡i21sin8⁡i-112sin6⁡i+240sin4⁡i-256sin2⁡i+128cos⁡2f+g-6365536sin2⁡i-263sin8⁡i-224sin6⁡i+352sin4⁡i-256sin2⁡i+128,(46)P10x=46189256x10-109395256x8+45045128x6-15015128x4+3465256x2-63256,(47)P10l′=290990733554432sin10⁡icos⁡10f+g-76576516777216sin8⁡i19sin2⁡i-18cos⁡8f+g+40540533554432sin6⁡i323sin4⁡i-544sin2⁡i+224cos⁡6f+g-1351354194304sin4⁡i323sin6⁡i-714sin4⁡i+504sin2⁡i-112cos⁡4f+g+7276516777216sin2⁡i4199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2f+g-18332414116777216sin10⁡i+2412159758388608sin8⁡i-14189175524288sin6⁡i+2837835262144sin4⁡i-218295131072sin2⁡i+396965536,(48)R10l′l=μ′a′aa′102565926893538589934592e10sin10⁡icos⁡10g-96463417054294967296e8sin8⁡ie2+619sin2⁡i-18cos⁡8g+895945058589934592e6sin6⁡i15e4+160e2+224323sin4⁡i-544sin2⁡i+224cos⁡6g-368918551073741824e4sin4⁡i5e6+70e4+168e2+80323sin6⁡i-714sin4⁡i+504sin2⁡i-112cos⁡4g+28378354294967296e2sin2⁡i21e8+336e6+1008e4+768e2+1284199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2g-39694294967296693e10+11550e8+36960e6+31680e4+7040e2+25646189sin10⁡i-121550sin8⁡i+114400sin6⁡i-45760sin4⁡i+7040sin2⁡i-256,(49)R10′l′l=μ′a′a′a111-e2-19/2-68918854294967296e8sin8⁡i19sin2⁡i-18cos⁡8g+12162151073741824e6sin6⁡i3e2+14323sin4⁡i-544sin2⁡i+224cos⁡6g-8513505268435456e4sin4⁡ie4+8e2+8323sin6i-714sin4⁡i+504sin2⁡i-112cos⁡4g+654885536870912e2sin2⁡i7e6+70e4+112e2+324199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2g-39692147483648315e8+3360e6+6048e4+2304e2+12846189sin10⁡i-121550sin8⁡i+114400sin6⁡i-45760sin4⁡i+7040sin2⁡i-256.


### 3.7.2n=12

At this order, the corresponding component of the disturbing function for the inner problemR12 is of O(α12), and that for the outer problem R12′ is of O(α′13).(50)
cos
12
⁡
S
l
′
=
231
2097152
sin
12
⁡
i
cos
⁡
12
f
+
g
-
693
524288
sin
10
⁡
i
sin
2
⁡
i
-
2
cos
⁡
10
f
+
g
+
693
1048576
sin
8
⁡
i
11
sin
4
⁡
i
-
40
sin
2
⁡
i
+
40
cos
⁡
8
f
+
g
-
1155
524288
sin
6
⁡
i
sin
2
⁡
i
-
2
11
sin
4
⁡
i
-
32
sin
2
⁡
i
+
32
cos
⁡
6
f
+
g
+
3465
2097152
sin
4
⁡
i
33
sin
8
⁡
i
-
192
sin
6
⁡
i
+
448
sin
4
⁡
i
-
512
sin
2
⁡
i
+
256
cos
⁡
4
f
+
g
-
693
262144
sin
2
⁡
i
sin
2
⁡
i
-
2
33
sin
8
⁡
i
-
144
sin
6
⁡
i
+
272
sin
4
⁡
i
-
256
sin
2
⁡
i
+
128
cos
⁡
2
f
+
g
+
53361
1048576
sin
12
⁡
i
-
43659
131072
sin
10
⁡
i
+
121275
131072
sin
8
⁡
i
-
5775
4096
sin
6
⁡
i
+
10395
8192
sin
4
⁡
i
-
693
1024
sin
2
⁡
i
+
231
1024
,
P
12
x
=
676039
1024
x
12
-
969969
512
x
10
+
2078505
1024
x
8
-
255255
256
x
6
+
225225
1024
x
4
-
9009
512
x
2
+
231
1024
,
P
12
l
′
=
156165009
2147483648
sin
12
⁡
i
cos
⁡
12
f
+
g
-
20369349
536870912
sin
10
⁡
i
23
sin
2
⁡
i
-
22
cos
⁡
10
f
+
g
+
32008977
1073741824
sin
8
⁡
i
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
f
+
g
-
2807805
536870912
sin
6
⁡
i
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
f
+
g
+
10405395
2147483648
sin
4
⁡
i
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
f
+
g
-
2081079
268435456
sin
2
⁡
i
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
f
+
g
+
36074117079
1073741824
sin
12
⁡
i
-
14115958857
134217728
sin
10
⁡
i
+
16804712925
134217728
sin
8
⁡
i
-
294819525
4194304
sin
6
⁡
i
+
156080925
8388608
sin
4
⁡
i
-
2081079
1048576
sin
2
⁡
i
+
53361
1048576
,
R
12
l
′
l
=
μ
′
a
′
a
a
′
12
203026224075675
2199023255552
e
12
sin
12
⁡
i
cos
⁡
12
g
-
1059267256047
549755813888
e
10
sin
10
⁡
i
3
e
2
+
22
23
sin
2
⁡
i
-
22
cos
⁡
10
g
+
361861484985
1099511627776
e
8
sin
8
⁡
i
3
e
4
+
40
e
2
+
72
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
g
-
6348447105
549755813888
e
6
sin
6
⁡
i
5
e
6
+
90
e
4
+
288
e
2
+
192
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
g
+
1238242005
2199023255552
e
4
sin
4
⁡
i
45
e
8
+
960
e
6
+
4032
e
4
+
4608
e
2
+
1280
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
g
-
72837765
274877906944
e
2
sin
2
⁡
i
99
e
10
+
2310
e
8
+
11088
e
6
+
15840
e
4
+
7040
e
2
+
768
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
g
+
53361
1099511627776
3003
e
12
+
72072
e
10
+
360360
e
8
+
549120
e
6
+
274560
e
4
+
39936
e
2
+
1024
676039
sin
12
⁡
i
-
2116296
sin
10
⁡
i
+
2519400
sin
8
⁡
i
-
1414400
sin
6
⁡
i
+
374400
sin
4
⁡
i
-
39936
sin
2
⁡
i
+
1024
,
R
12
′
l
′
l
=
μ
′
a
′
a
′
a
13
1
-
e
2
-
23
/
2
-
224062839
549755813888
e
10
sin
10
⁡
i
23
sin
2
⁡
i
-
22
cos
⁡
10
g
+
1760493735
549755813888
e
8
sin
8
⁡
i
e
2
+
6
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
g
-
92657565
549755813888
e
6
sin
6
⁡
i
15
e
4
+
160
e
2
+
224
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
g
+
343378035
274877906944
e
4
sin
4
⁡
i
5
e
6
+
70
e
4
+
168
e
2
+
80
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
g
-
114459345
137438953472
e
2
sin
2
⁡
i
21
e
8
+
336
e
6
+
1008
e
4
+
768
e
2
+
128
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
g
+
53361
274877906944
693
e
10
+
11550
e
8
+
36960
e
6
+
31680
e
4
+
7040
e
2
+
256
676039
sin
12
⁡
i
-
2116296
sin
10
⁡
i
+
2519400
sin
8
⁡
i
-
1414400
sin
6
⁡
i
+
374400
sin
4
⁡
i
-
39936
sin
2
⁡
i
+
1024
.


### 3.8.2n=14

At this order, the corresponding component of the disturbing function for the inner problemR14 is of O(α14), and that for the outer problem R14′ is of O(α′15).(51)
cos
14
⁡
S
l
′
=
429
16777216
sin
14
⁡
i
cos
⁡
14
f
+
g
-
3003
8388608
sin
12
⁡
i
sin
2
⁡
i
-
2
cos
⁡
12
f
+
g
+
3003
16777216
sin
10
⁡
i
13
sin
4
⁡
i
-
48
sin
2
⁡
i
+
48
cos
⁡
10
f
+
g
-
3003
4194304
sin
8
⁡
i
sin
2
⁡
i
-
2
13
sin
4
⁡
i
-
40
sin
2
⁡
i
+
40
cos
⁡
8
f
+
g
+
3003
16777216
sin
6
⁡
i
143
sin
8
⁡
i
-
880
sin
6
⁡
i
+
2160
sin
4
⁡
i
-
2560
sin
2
⁡
i
+
1280
cos
⁡
6
f
+
g
-
3003
8388608
sin
4
⁡
i
sin
2
⁡
i
-
2
143
sin
8
⁡
i
-
704
sin
6
⁡
i
+
1472
sin
4
⁡
i
-
1536
sin
2
⁡
i
+
768
cos
⁡
4
f
+
g
+
3003
16777216
sin
2
⁡
i
429
sin
12
⁡
i
-
3168
sin
10
⁡
i
+
10080
sin
8
⁡
i
-
17920
sin
6
⁡
i
+
19200
sin
4
⁡
i
-
12288
sin
2
⁡
i
+
4096
cos
⁡
2
f
+
g
-
429
4194304
sin
2
⁡
i
-
2
429
sin
12
⁡
i
-
2376
sin
10
⁡
i
+
5832
sin
8
⁡
i
-
7936
sin
6
⁡
i
+
6528
sin
4
⁡
i
-
3072
sin
2
⁡
i
+
1024
,
P
14
x
=
5014575
2048
x
14
-
16900975
2048
x
12
+
22309287
2048
x
10
-
14549535
2048
x
8
+
4849845
2048
x
6
-
765765
2048
x
4
+
45045
2048
x
2
-
429
2048
,
P
14
l
′
=
2151252675
34359738368
sin
14
⁡
i
cos
⁡
14
f
+
g
-
557732175
17179869184
sin
12
⁡
i
27
sin
2
⁡
i
-
26
cos
⁡
12
f
+
g
+
870062193
34359738368
sin
10
⁡
i
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
f
+
g
-
189143955
8589934592
sin
8
⁡
i
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
f
+
g
+
693527835
34359738368
sin
6
⁡
i
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
f
+
g
-
328513185
17179869184
sin
4
⁡
i
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
f
+
g
+
19324305
34359738368
sin
2
⁡
i
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
f
+
g
-
922887397575
8589934592
sin
14
⁡
i
+
1674869721525
4294967296
sin
12
⁡
i
-
602953099749
1073741824
sin
10
⁡
i
+
218461268025
536870912
sin
8
⁡
i
-
10402917525
67108864
sin
6
⁡
i
+
985539555
33554432
sin
4
⁡
i
-
19324305
8388608
sin
2
⁡
i
+
184041
4194304
,
R
14
l
′
l
=
μ
′
a
′
a
a
′
14
20856061239960375
70368744177664
e
14
sin
14
⁡
i
cos
⁡
14
g
-
186452654763375
35184372088832
e
12
sin
12
⁡
i
3
e
2
+
26
27
sin
2
⁡
i
-
26
cos
⁡
12
g
+
32318460158985
70368744177664
e
10
sin
10
⁡
i
5
e
4
+
80
e
2
+
176
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
g
-
1405150441695
17592186044416
e
8
sin
8
⁡
i
7
e
6
+
154
e
4
+
616
e
2
+
528
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
g
+
672028472115
70368744177664
e
6
sin
6
⁡
i
21
e
8
+
560
e
6
+
3024
e
4
+
4608
e
2
+
1792
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
g
-
318329276265
35184372088832
e
4
sin
4
⁡
i
11
e
10
+
330
e
8
+
2112
e
6
+
4224
e
4
+
2816
e
2
+
512
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
g
+
328513185
70368744177664
e
2
sin
2
⁡
i
429
e
12
+
13728
e
10
+
96096
e
8
+
219648
e
6
+
183040
e
4
+
53248
e
2
+
4096
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
g
-
184041
17592186044416
6435
e
14
+
210210
e
12
+
1513512
e
10
+
3603600
e
8
+
3203200
e
6
+
1048320
e
4
+
107520
e
2
+
2048
5014575
sin
14
⁡
i
-
18201050
sin
12
⁡
i
+
26209512
sin
10
⁡
i
-
18992400
sin
8
⁡
i
+
7235200
sin
6
⁡
i
-
1370880
sin
4
⁡
i
+
107520
sin
2
⁡
i
-
2048
,
R
14
′
l
′
l
=
μ
′
a
′
a
′
a
15
1
-
e
2
-
27
/
2
-
7250518275
70368744177664
e
12
sin
12
⁡
i
27
sin
2
⁡
i
-
26
cos
⁡
12
g
+
11310808509
35184372088832
e
10
sin
10
⁡
i
3
e
2
+
22
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
g
-
27047585565
17592186044416
e
8
sin
8
⁡
i
3
e
4
+
40
e
2
+
72
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
g
+
99174480405
35184372088832
e
6
sin
6
⁡
i
5
e
6
+
90
e
4
+
288
e
2
+
192
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
g
-
46977385455
70368744177664
e
4
sin
4
⁡
i
45
e
8
+
960
e
6
+
4032
e
4
+
4608
e
2
+
1280
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
g
+
251215965
17592186044416
e
2
sin
2
⁡
i
99
e
10
+
2310
e
8
+
11088
e
6
+
15840
e
4
+
7040
e
2
+
768
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
g
-
184041
8796093022208
3003
e
12
+
72072
e
10
+
360360
e
8
+
549120
e
6
+
274560
e
4
+
39936
e
2
+
1024
5014575
sin
14
⁡
i
-
18201050
sin
12
⁡
i
+
26209512
sin
10
⁡
i
-
18992400
sin
8
⁡
i
+
7235200
sin
6
⁡
i
-
1370880
sin
4
⁡
i
+
107520
sin
2
⁡
i
-
2048
.


## 3.1. General Form

From (3) and (4), the direct part of the disturbing function for the inner CR3BP where r′=a′ becomes as follows:(12)R=μ′a′∑n=1∞ra′2nP2ncos⁡S,where we ignore the term of n=0. In (12) we also ignore all the terms including the odd Legendre polynomials (P1,P3,P5,P7,P9,…) 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 R 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 ofR (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 l and l′ in what follows) are independent from each other. The procedure to carry out double averaging of R is straightforward as follows: Let us pick up the nth term of R and name it as R2n. We have(13)R2n=μ′a′ra′2nP2ncos⁡S.First we averageR2n by mean anomaly of the perturbing body l′, as(14)R2nl′=μ′a′ra′2nP2nl′,where(15)P2nl′=12π∫02πP2ncos⁡Sdl′.The angleS is expressed by orbital angles through a relationship (e.g., [3, Eq. (7) in p. 592]) (16)cos⁡S=cos⁡f+gcos⁡f′+g′+cos⁡isin⁡f+gsin⁡f′+g′,where f, f′ are true anomalies of the perturbed and perturbing bodies, g, g′ are arguments of pericenter of the perturbed and perturbing bodies, and i 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 g and g′ can be measured from the mutual node. Note that g′ is not actually defined in CR3BP. Therefore, in (16) we regard f′+g′ as a single, fast-moving variable when we carry out averaging of (15). Practically, we can simply replace ∫dl′ for ∫df′ in the discussion here.To obtainP2nl′ of (15), we calculate the average of cos2n⁡S by l′ as(17)cos2n⁡Sl′=12π∫02πcos2n⁡Sdl′.Then we average R2nl′ of (14) by mean anomaly of the perturbed body l, as(18)R2nl′l=μ′a′aa′2n12π∫02πra2nP2nl′dl.If we switch the integration variable froml to eccentric anomaly u, (18) becomes(19)R2nl′l=μ′a′aa′2n12π∫02π1-ecos⁡u2n+1P2nl′du.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 R′ for the outer CR3BP becomes as follows:(20)R′=μ′r∑n=1∞a′r2nP2ncos⁡S.Note that our definition of1/Δ for the outer case (5), hence also in (20), may be different from what is seen in conventional textbooks (e.g., [15, Eq. (6.22) in p. 229]): Roles of the dashed quantities (r′,μ′) may be the opposite. This difference comes from the fact that conventional textbooks always assume r/r′<1, while we assume r/r′>1 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 ofR′ over mean anomalies of both the bodies. Let us pick the nth term of R′ in (20) and name it R2n′. We have(21)R2n′=μ′ra′r2nP2ncos⁡S.First we averageR2n′ by mean anomaly of the perturbing body l′, as(22)R2n′l′=μ′ra′r2nP2nl′,where P2nl′ is already defined in (15).Then we averageR2n′l′ in (22) by mean anomaly of the perturbed body l, as(23)R2n′l′l=μ′a′a′a2n+112π∫02πar2n+1P2nl′dl.If we switch the integration variable froml to true anomaly f, (23) becomes(24)R2n′l′l=μ′a′a′a2n+11-e2-2n+1/22π∫02π1+ecos⁡f2n-1P2nl′df.Note thatR2n′l′l has the order of O(α′2n+1) as in (23) and (24), not O(α′2n).In Sections3.2 to 3.8 we show the expressions of cos2n⁡Sl′ in (17), P2nl′ in (15), the Legendre polynomial P2n(x) in its original form, R2nl′l in (19), and R2n′l′l in (24) for 2n=2,4,6,8,10,12,14. We used Maple™ for algebraic manipulation to obtain the series of expressions. Note that in what follows we use sin⁡i instead of cos⁡i 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.2n=2

At this order, the corresponding component of the disturbing function for the inner problemR2 is of O(α2), and that for the outer problem R2′ is of O(α′3).We just describe the resulting expressions of the expansion as follows: Let us emphasize again that the dashed quantities such asμ′ and a′ are those of the perturbing body, whether its orbit is located inside or outside of the orbit of the perturbed body.(25)cos2⁡Sl′=14sin2⁡icos⁡2f+g-14sin2⁡i+12,(26)P2x=32x2-12,(27)P2l′=38sin2⁡icos⁡2f+g-38sin2⁡i+14,(28)R2l′l=μ′a′aa′21516e2sin2⁡icos⁡2g-1163e2+23sin2⁡i-2,(29)R2′l′l=μ′a′a′a31-e2-3/2-38sin2⁡i+14.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 g. The g-dependence in the doubly averaged outer CR3BP first shows up in the next order: 2n=4. 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.2n=4

At this order, the corresponding component of the disturbing function for the inner problemR4 is of O(α4), and that for the outer problem R4′ is of O(α′5).(30)cos4⁡Sl′=364sin4⁡icos⁡4f+g-316sin2⁡isin2⁡i-2cos⁡2f+g+964sin4⁡i-38sin2⁡i+38,(31)P4x=358x4-154x2+38,(32)P4l′=105512sin4⁡icos⁡4f+g-15128sin2⁡i7sin2⁡i-6cos⁡2f+g+315512sin4⁡i-4564sin2⁡i+964,(33)R4l′l=μ′a′aa′466154096e4sin4⁡icos⁡4g-3151024e2sin2⁡ie2+27sin2⁡i-6cos⁡2g+9409615e4+40e2+835sin4⁡i+40sin2⁡i+8,(34)R4′l′l=μ′a′a′a51-e2-7/2-45512e2sin2⁡i7sin2⁡i-6cos⁡2g+910243e2+235sin4⁡i-40sin2⁡i+8.Note that we now see theg-dependence of the doubly averaged disturbing function for the outer CR3BP in the expression of (34).


## 3.4.2n=6

At this order, the corresponding component of the disturbing function for the inner problemR6 is of O(α6), and that for the outer problem R6′ is of O(α′7).(35)cos6⁡Sl′=5512sin6⁡icos⁡6f+g-15256sin4⁡isin2⁡i-2cos⁡4f+g+15512sin2⁡i5sin4⁡i-16sin2⁡i+16cos⁡2f+g-5256sin2⁡i-25sin4⁡i-8sin2⁡i+8,(36)P6x=23116x6-31516x4+10516x2-516,(37)P6l′=11558192sin6⁡icos⁡6f+g-3154096sin4⁡i11sin2⁡i-10cos⁡4f+g+5258192sin2⁡i33sin4⁡i-48sin2⁡i+16cos⁡2f+g-57754096sin6⁡i+47252048sin4⁡i-525512sin2⁡i+25256,(38)R6l′l=μ′a′aa′6495495131072e6sin6⁡icos⁡6g-1039565536e4sin4⁡i3e2+1011sin2⁡i-10cos⁡4g+1575131072e2sin2⁡i15e4+80e2+4833sin4⁡i-48sin2⁡i+16cos⁡2g-256553635e6+210e4+168e2+16231sin6⁡i-378sin4⁡i+168sin2⁡i-16,(39)R6′l′l=μ′a′a′a71-e2-11/2-157565536e4sin4⁡i11sin2⁡i-10cos⁡4g+262532768e2sin2⁡ie2+233sin4⁡i-48sin2⁡i+16cos⁡2g-253276815e4+40e2+8231sin6⁡i-378sin4⁡i+168sin2⁡i-16.


## 3.5.2n=8

At this order, the corresponding component of the disturbing function for the inner problemR8 is of O(α8), and that for the outer problem R8′ is of O(α′9).(40)cos8⁡Sl′=3516384sin8⁡icos⁡8f+g-352048sin6⁡isin2⁡i-2cos⁡6f+g+354096sin4⁡i7sin4⁡i-24sin2⁡i+24cos⁡4f+g-352048sin2⁡isin2⁡i-27sin4⁡i-16sin2⁡i+16cos⁡2f+g+122516384sin8⁡i-175512sin6⁡i+315512sin4⁡i-3564sin2⁡i+35128,(41)P8x=6435128x8-300332x6+346564x4-31532x2+35128,(42)P8l′=2252252097152sin8⁡icos⁡8f+g-15015262144sin6⁡i15sin2⁡i-14cos⁡6f+g+24255524288sin4⁡i65sin4⁡i-104sin2⁡i+40cos⁡4f+g-11025262144sin2⁡i143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2f+g+78828752097152sin8⁡i-52552565536sin6⁡i+36382565536sin4⁡i-110258192sin2⁡i+122516384,(43)R8l′l=μ′a′aa′82737609875268435456e8sin8⁡icos⁡8g-1073572533554432e6sin6⁡i3e2+1415sin2⁡i-14cos⁡6g+1734232567108864e4sin4⁡ie4+8e2+865sin4⁡i-104sin2⁡i+40cos⁡4g-60637533554432e2sin2⁡i7e6+70e4+112e2+32143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2g+1225268435456315e8+3360e6+6048e4+2304e2+1286435sin8⁡i-13728sin6⁡i+9504sin4⁡i-2304sin2⁡i+128,(44)R8′l′l=μ′a′a′a91-e2-15/2-10510516777216e6sin6⁡i15sin2⁡i-14cos⁡6g+16978516777216e4sin4⁡i3e2+1065sin4⁡i-104sin2⁡i+40cos⁡4g-7717516777216e2sin2⁡i15e4+80e2+48143sin6⁡i-286sin4⁡i+176sin2⁡i-32cos⁡2g+12253355443235e6+210e4+168e2+166435sin8⁡i-13728sin6⁡i+9504sin4⁡i-2304sin2⁡i+128.


## 3.6.2n=10

At this order, the corresponding component of the disturbing function for the inner problemR10 is of O(α10), and that for the outer problem R10′ is of O(α′11).(45)cos10⁡Sl′=63131072sin10⁡icos⁡10f+g-31565536sin8⁡isin2⁡i-2cos⁡8f+g+315131072sin6⁡i9sin4⁡i-32sin2⁡i+32cos⁡6f+g-31516384sin4⁡isin2⁡i-23sin4⁡i-8sin2⁡i+8cos⁡4f+g+31565536sin2⁡i21sin8⁡i-112sin6⁡i+240sin4⁡i-256sin2⁡i+128cos⁡2f+g-6365536sin2⁡i-263sin8⁡i-224sin6⁡i+352sin4⁡i-256sin2⁡i+128,(46)P10x=46189256x10-109395256x8+45045128x6-15015128x4+3465256x2-63256,(47)P10l′=290990733554432sin10⁡icos⁡10f+g-76576516777216sin8⁡i19sin2⁡i-18cos⁡8f+g+40540533554432sin6⁡i323sin4⁡i-544sin2⁡i+224cos⁡6f+g-1351354194304sin4⁡i323sin6⁡i-714sin4⁡i+504sin2⁡i-112cos⁡4f+g+7276516777216sin2⁡i4199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2f+g-18332414116777216sin10⁡i+2412159758388608sin8⁡i-14189175524288sin6⁡i+2837835262144sin4⁡i-218295131072sin2⁡i+396965536,(48)R10l′l=μ′a′aa′102565926893538589934592e10sin10⁡icos⁡10g-96463417054294967296e8sin8⁡ie2+619sin2⁡i-18cos⁡8g+895945058589934592e6sin6⁡i15e4+160e2+224323sin4⁡i-544sin2⁡i+224cos⁡6g-368918551073741824e4sin4⁡i5e6+70e4+168e2+80323sin6⁡i-714sin4⁡i+504sin2⁡i-112cos⁡4g+28378354294967296e2sin2⁡i21e8+336e6+1008e4+768e2+1284199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2g-39694294967296693e10+11550e8+36960e6+31680e4+7040e2+25646189sin10⁡i-121550sin8⁡i+114400sin6⁡i-45760sin4⁡i+7040sin2⁡i-256,(49)R10′l′l=μ′a′a′a111-e2-19/2-68918854294967296e8sin8⁡i19sin2⁡i-18cos⁡8g+12162151073741824e6sin6⁡i3e2+14323sin4⁡i-544sin2⁡i+224cos⁡6g-8513505268435456e4sin4⁡ie4+8e2+8323sin6i-714sin4⁡i+504sin2⁡i-112cos⁡4g+654885536870912e2sin2⁡i7e6+70e4+112e2+324199sin8⁡i-10608sin6⁡i+9360sin4⁡i-3328sin2⁡i+384cos⁡2g-39692147483648315e8+3360e6+6048e4+2304e2+12846189sin10⁡i-121550sin8⁡i+114400sin6⁡i-45760sin4⁡i+7040sin2⁡i-256.


## 3.7.2n=12

At this order, the corresponding component of the disturbing function for the inner problemR12 is of O(α12), and that for the outer problem R12′ is of O(α′13).(50)
cos
12
⁡
S
l
′
=
231
2097152
sin
12
⁡
i
cos
⁡
12
f
+
g
-
693
524288
sin
10
⁡
i
sin
2
⁡
i
-
2
cos
⁡
10
f
+
g
+
693
1048576
sin
8
⁡
i
11
sin
4
⁡
i
-
40
sin
2
⁡
i
+
40
cos
⁡
8
f
+
g
-
1155
524288
sin
6
⁡
i
sin
2
⁡
i
-
2
11
sin
4
⁡
i
-
32
sin
2
⁡
i
+
32
cos
⁡
6
f
+
g
+
3465
2097152
sin
4
⁡
i
33
sin
8
⁡
i
-
192
sin
6
⁡
i
+
448
sin
4
⁡
i
-
512
sin
2
⁡
i
+
256
cos
⁡
4
f
+
g
-
693
262144
sin
2
⁡
i
sin
2
⁡
i
-
2
33
sin
8
⁡
i
-
144
sin
6
⁡
i
+
272
sin
4
⁡
i
-
256
sin
2
⁡
i
+
128
cos
⁡
2
f
+
g
+
53361
1048576
sin
12
⁡
i
-
43659
131072
sin
10
⁡
i
+
121275
131072
sin
8
⁡
i
-
5775
4096
sin
6
⁡
i
+
10395
8192
sin
4
⁡
i
-
693
1024
sin
2
⁡
i
+
231
1024
,
P
12
x
=
676039
1024
x
12
-
969969
512
x
10
+
2078505
1024
x
8
-
255255
256
x
6
+
225225
1024
x
4
-
9009
512
x
2
+
231
1024
,
P
12
l
′
=
156165009
2147483648
sin
12
⁡
i
cos
⁡
12
f
+
g
-
20369349
536870912
sin
10
⁡
i
23
sin
2
⁡
i
-
22
cos
⁡
10
f
+
g
+
32008977
1073741824
sin
8
⁡
i
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
f
+
g
-
2807805
536870912
sin
6
⁡
i
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
f
+
g
+
10405395
2147483648
sin
4
⁡
i
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
f
+
g
-
2081079
268435456
sin
2
⁡
i
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
f
+
g
+
36074117079
1073741824
sin
12
⁡
i
-
14115958857
134217728
sin
10
⁡
i
+
16804712925
134217728
sin
8
⁡
i
-
294819525
4194304
sin
6
⁡
i
+
156080925
8388608
sin
4
⁡
i
-
2081079
1048576
sin
2
⁡
i
+
53361
1048576
,
R
12
l
′
l
=
μ
′
a
′
a
a
′
12
203026224075675
2199023255552
e
12
sin
12
⁡
i
cos
⁡
12
g
-
1059267256047
549755813888
e
10
sin
10
⁡
i
3
e
2
+
22
23
sin
2
⁡
i
-
22
cos
⁡
10
g
+
361861484985
1099511627776
e
8
sin
8
⁡
i
3
e
4
+
40
e
2
+
72
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
g
-
6348447105
549755813888
e
6
sin
6
⁡
i
5
e
6
+
90
e
4
+
288
e
2
+
192
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
g
+
1238242005
2199023255552
e
4
sin
4
⁡
i
45
e
8
+
960
e
6
+
4032
e
4
+
4608
e
2
+
1280
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
g
-
72837765
274877906944
e
2
sin
2
⁡
i
99
e
10
+
2310
e
8
+
11088
e
6
+
15840
e
4
+
7040
e
2
+
768
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
g
+
53361
1099511627776
3003
e
12
+
72072
e
10
+
360360
e
8
+
549120
e
6
+
274560
e
4
+
39936
e
2
+
1024
676039
sin
12
⁡
i
-
2116296
sin
10
⁡
i
+
2519400
sin
8
⁡
i
-
1414400
sin
6
⁡
i
+
374400
sin
4
⁡
i
-
39936
sin
2
⁡
i
+
1024
,
R
12
′
l
′
l
=
μ
′
a
′
a
′
a
13
1
-
e
2
-
23
/
2
-
224062839
549755813888
e
10
sin
10
⁡
i
23
sin
2
⁡
i
-
22
cos
⁡
10
g
+
1760493735
549755813888
e
8
sin
8
⁡
i
e
2
+
6
161
sin
4
⁡
i
-
280
sin
2
⁡
i
+
120
cos
⁡
8
g
-
92657565
549755813888
e
6
sin
6
⁡
i
15
e
4
+
160
e
2
+
224
3059
sin
6
⁡
i
-
7182
sin
4
⁡
i
+
5472
sin
2
⁡
i
-
1344
cos
⁡
6
g
+
343378035
274877906944
e
4
sin
4
⁡
i
5
e
6
+
70
e
4
+
168
e
2
+
80
7429
sin
8
⁡
i
-
20672
sin
6
⁡
i
+
20672
sin
4
⁡
i
-
8704
sin
2
⁡
i
+
1280
cos
⁡
4
g
-
114459345
137438953472
e
2
sin
2
⁡
i
21
e
8
+
336
e
6
+
1008
e
4
+
768
e
2
+
128
7429
sin
10
⁡
i
-
22610
sin
8
⁡
i
+
25840
sin
6
⁡
i
-
13600
sin
4
⁡
i
+
3200
sin
2
⁡
i
-
256
cos
⁡
2
g
+
53361
274877906944
693
e
10
+
11550
e
8
+
36960
e
6
+
31680
e
4
+
7040
e
2
+
256
676039
sin
12
⁡
i
-
2116296
sin
10
⁡
i
+
2519400
sin
8
⁡
i
-
1414400
sin
6
⁡
i
+
374400
sin
4
⁡
i
-
39936
sin
2
⁡
i
+
1024
.


## 3.8.2n=14

At this order, the corresponding component of the disturbing function for the inner problemR14 is of O(α14), and that for the outer problem R14′ is of O(α′15).(51)
cos
14
⁡
S
l
′
=
429
16777216
sin
14
⁡
i
cos
⁡
14
f
+
g
-
3003
8388608
sin
12
⁡
i
sin
2
⁡
i
-
2
cos
⁡
12
f
+
g
+
3003
16777216
sin
10
⁡
i
13
sin
4
⁡
i
-
48
sin
2
⁡
i
+
48
cos
⁡
10
f
+
g
-
3003
4194304
sin
8
⁡
i
sin
2
⁡
i
-
2
13
sin
4
⁡
i
-
40
sin
2
⁡
i
+
40
cos
⁡
8
f
+
g
+
3003
16777216
sin
6
⁡
i
143
sin
8
⁡
i
-
880
sin
6
⁡
i
+
2160
sin
4
⁡
i
-
2560
sin
2
⁡
i
+
1280
cos
⁡
6
f
+
g
-
3003
8388608
sin
4
⁡
i
sin
2
⁡
i
-
2
143
sin
8
⁡
i
-
704
sin
6
⁡
i
+
1472
sin
4
⁡
i
-
1536
sin
2
⁡
i
+
768
cos
⁡
4
f
+
g
+
3003
16777216
sin
2
⁡
i
429
sin
12
⁡
i
-
3168
sin
10
⁡
i
+
10080
sin
8
⁡
i
-
17920
sin
6
⁡
i
+
19200
sin
4
⁡
i
-
12288
sin
2
⁡
i
+
4096
cos
⁡
2
f
+
g
-
429
4194304
sin
2
⁡
i
-
2
429
sin
12
⁡
i
-
2376
sin
10
⁡
i
+
5832
sin
8
⁡
i
-
7936
sin
6
⁡
i
+
6528
sin
4
⁡
i
-
3072
sin
2
⁡
i
+
1024
,
P
14
x
=
5014575
2048
x
14
-
16900975
2048
x
12
+
22309287
2048
x
10
-
14549535
2048
x
8
+
4849845
2048
x
6
-
765765
2048
x
4
+
45045
2048
x
2
-
429
2048
,
P
14
l
′
=
2151252675
34359738368
sin
14
⁡
i
cos
⁡
14
f
+
g
-
557732175
17179869184
sin
12
⁡
i
27
sin
2
⁡
i
-
26
cos
⁡
12
f
+
g
+
870062193
34359738368
sin
10
⁡
i
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
f
+
g
-
189143955
8589934592
sin
8
⁡
i
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
f
+
g
+
693527835
34359738368
sin
6
⁡
i
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
f
+
g
-
328513185
17179869184
sin
4
⁡
i
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
f
+
g
+
19324305
34359738368
sin
2
⁡
i
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
f
+
g
-
922887397575
8589934592
sin
14
⁡
i
+
1674869721525
4294967296
sin
12
⁡
i
-
602953099749
1073741824
sin
10
⁡
i
+
218461268025
536870912
sin
8
⁡
i
-
10402917525
67108864
sin
6
⁡
i
+
985539555
33554432
sin
4
⁡
i
-
19324305
8388608
sin
2
⁡
i
+
184041
4194304
,
R
14
l
′
l
=
μ
′
a
′
a
a
′
14
20856061239960375
70368744177664
e
14
sin
14
⁡
i
cos
⁡
14
g
-
186452654763375
35184372088832
e
12
sin
12
⁡
i
3
e
2
+
26
27
sin
2
⁡
i
-
26
cos
⁡
12
g
+
32318460158985
70368744177664
e
10
sin
10
⁡
i
5
e
4
+
80
e
2
+
176
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
g
-
1405150441695
17592186044416
e
8
sin
8
⁡
i
7
e
6
+
154
e
4
+
616
e
2
+
528
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
g
+
672028472115
70368744177664
e
6
sin
6
⁡
i
21
e
8
+
560
e
6
+
3024
e
4
+
4608
e
2
+
1792
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
g
-
318329276265
35184372088832
e
4
sin
4
⁡
i
11
e
10
+
330
e
8
+
2112
e
6
+
4224
e
4
+
2816
e
2
+
512
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
g
+
328513185
70368744177664
e
2
sin
2
⁡
i
429
e
12
+
13728
e
10
+
96096
e
8
+
219648
e
6
+
183040
e
4
+
53248
e
2
+
4096
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
g
-
184041
17592186044416
6435
e
14
+
210210
e
12
+
1513512
e
10
+
3603600
e
8
+
3203200
e
6
+
1048320
e
4
+
107520
e
2
+
2048
5014575
sin
14
⁡
i
-
18201050
sin
12
⁡
i
+
26209512
sin
10
⁡
i
-
18992400
sin
8
⁡
i
+
7235200
sin
6
⁡
i
-
1370880
sin
4
⁡
i
+
107520
sin
2
⁡
i
-
2048
,
R
14
′
l
′
l
=
μ
′
a
′
a
′
a
15
1
-
e
2
-
27
/
2
-
7250518275
70368744177664
e
12
sin
12
⁡
i
27
sin
2
⁡
i
-
26
cos
⁡
12
g
+
11310808509
35184372088832
e
10
sin
10
⁡
i
3
e
2
+
22
45
sin
2
⁡
i
-
44
5
sin
2
⁡
i
-
4
cos
⁡
10
g
-
27047585565
17592186044416
e
8
sin
8
⁡
i
3
e
4
+
40
e
2
+
72
1035
sin
6
⁡
i
-
2530
sin
4
⁡
i
+
2024
sin
2
⁡
i
-
528
cos
⁡
8
g
+
99174480405
35184372088832
e
6
sin
6
⁡
i
5
e
6
+
90
e
4
+
288
e
2
+
192
3105
sin
8
⁡
i
-
9200
sin
6
⁡
i
+
9936
sin
4
⁡
i
-
4608
sin
2
⁡
i
+
768
cos
⁡
6
g
-
46977385455
70368744177664
e
4
sin
4
⁡
i
45
e
8
+
960
e
6
+
4032
e
4
+
4608
e
2
+
1280
6555
sin
10
⁡
i
-
21850
sin
8
⁡
i
+
27968
sin
6
⁡
i
-
17024
sin
4
⁡
i
+
4864
sin
2
⁡
i
-
512
cos
⁡
4
g
+
251215965
17592186044416
e
2
sin
2
⁡
i
99
e
10
+
2310
e
8
+
11088
e
6
+
15840
e
4
+
7040
e
2
+
768
334305
sin
12
⁡
i
-
1188640
sin
10
⁡
i
+
1664096
sin
8
⁡
i
-
1157632
sin
6
⁡
i
+
413440
sin
4
⁡
i
-
69632
sin
2
⁡
i
+
4096
cos
⁡
2
g
-
184041
8796093022208
3003
e
12
+
72072
e
10
+
360360
e
8
+
549120
e
6
+
274560
e
4
+
39936
e
2
+
1024
5014575
sin
14
⁡
i
-
18201050
sin
12
⁡
i
+
26209512
sin
10
⁡
i
-
18992400
sin
8
⁡
i
+
7235200
sin
6
⁡
i
-
1370880
sin
4
⁡
i
+
107520
sin
2
⁡
i
-
2048
.


## 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)(52)Rl′l=μ′4π2∬02πdldl′Δ.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 (x,y)=ecos⁡g,esin⁡g plane. Since g 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 k2 defined as (e.g., [2, 34–37])(53)k2=1-e2cos2⁡i.We selected several combinations of(α,k2) 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 e is achieved when cos⁡i=1 as emax=1-k2. Similarly, the theoretical smallest value of cos⁡i is achieved when e=0 as cos⁡imin=k.Figure 2
Equipotential curves computed through the doubly averaged disturbing function and plotted on the plane(x,y)=(ecos⁡g,esin⁡g). Columns are, from the left, 2n=2 (for (a), (b)) or 2n=4 (for (c), (d), (e), and (f)), 2n=8, 2n=14, and the numerical quadrature defined by (52). Rows are, from the top, (a) (α,k2)=(0.90,0.900), (b) (α,k2)=(0.90,0.400), (c) (α′,k2)=(0.70,0.600), (d) (α′,k2)=(0.60,0.265), (e) (α′,k2)=(0.60,0.260), and (f) (α′,k2)=(0.80,0.304) for the analytic expansion and (α′,k2)=(0.80,0.313) for the numerical quadrature. The black dashed circles represent the theoretical maximum eccentricity (emax) of each of the systems. The red and the blue partial circles represent the conditions where the orbits of perturbed and perturbing bodies cross at the ascending node (red) and at the descending node (blue) of the perturbed body.
(a) (b) (c) (d) (e)  (f)Figures2(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 a and a′ (i.e., functions of α), and the relationship between a, e, g, and a′ is expressed as (e.g., [37–40])(54)a′=a1-e21±ecos⁡g.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 f=-g, where f denotes perturbed body’s true anomaly). The negative sign takes place when the encounter happens at the descending node of perturbed body (when f=π-g). Using the variables x=ecos⁡g and y=esin⁡g, (54) can be rewritten as(55)x±12α2+y2=1-12α2,and it is more obvious that (55) represents a pair of circles on the (x,y) plane.Note that the orbit-crossing lines introduce singularities inN-body Hamiltonian (e.g., [41, 42]), and R (or R′) 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 2n=14 and α=0.50, α14=6.103515625×10-5 which we may ignore. But when 2n=14 and α=0.90, α14=0.22876792454961, which is not ignorable. In Figure 2(a) where (α,k2)=(0.90,0.900), effect of the high-order analytic expansion is well exhibited. As is widely known, there is no stationary point for argument of pericenter g while k2>0.6 in the doubly averaged inner CR3BP at the quadruple-order approximation (2n=2). However when α becomes large, g acquires stationary points at g=±π/2 even when k2>0.6. 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 2n=8. But the stationary points of g seen in the rightmost panel of Figure 2(a) (for quadrature) are not reproduced even by using the 2n=8 analytic terms. We confirmed that the stationary points at g=±π/2 in this parameter set (α,k2)=(0.90,0.900) first appear in the analytic approximation up to 2n=12, and the value of e at the stationary points is better reproduced using the analytic expansion of 2n=14, as seen in Figure 2(a). Note that other sets of stationary points at g=0,π 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(e=emax) 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.Figure2(b) is for (α,k2)=(0.90,0.400). 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 (x,y)=(0,0) is reproduced in the results of the analytic expansion of disturbing function, we see spurious local extremums along the y-axis near the boundary of e=emax in the 2n=14 case. The spurious local extremums are not so remarkable in the 2n=8 case (although they surely exist), but there seem to exist spurious saddle points around y=esin⁡g~±0.3 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 2n=2 looks the closest to the quadrature result. However, we think this is just a coincidence, because the 2n=2 expansion includes only the terms of cos⁡2g as in (28), and R2l′l would not be influenced by terms of cos⁡2pg with p≥2.Figure2(c) is for the doubly averaged outer CR3BP when (α′,k2)=(0.70,0.600). In this case, the high-order analytic expansion of the doubly averaged disturbing function (2n=14) reproduces the quadrature result inside the orbit-crossing lines. However, the medium-order expansion (2n=8) yields artificial saddle points along the y-axis. We confirmed that these spurious saddle points show up in the expansion of 2n=10 too (panels are not shown here). We thus say that the highest-order analytic expansion of 2n=14 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 e is moderate.In Figure2(d) when (α′,k2)=(0.60,0.265), we see three local extremums in the quadrature result: a pair of local minima along the y-axis together with a saddle point at the origin (0,0). These local extremums cannot be reproduced by the low-order analytic expansion such as 2n=4 or 2n=8; the local extremums show up only in the expansion of 2n=12 or higher, although the locations of the local minima along the y-axis are slightly different from the quadrature result. When k2 becomes even smaller, the difference in the topology of the equipotential curves increases, as seen in Figure 2(e) where (α′,k2)=(0.60,0.260). A pair of local minima along the y-axis persists, but the origin (0,0) becomes a local maximum and a pair of saddle points takes place along the x-axis. We confirmed that this topology is only reproduced when using the analytic expansion of disturbing function at our highest-order (2n=14), although the locations of the local extremums are slightly different from the quadrature result.Whenα′ gets even larger and k2 is in a certain range of moderate values, a pair of local minima show up along the x-axis inside the orbit-crossing lines together with a saddle point at the origin (0,0). This is the case of Figure 2(f) when (α′,k2)=(0.80,0.304). The topology of the local extremums is reproduced by the analytic expansion of disturbing function at our highest-order (2n=14). 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 x-axis appear when the value of k2 is slightly different (k2=0.313) from the analytic result (k2=0.304). We have not yet figured out the cause of the discrepancy of the k2 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 (2n=16,18,20,…) 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 Section4 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 k2 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 Figure2, 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.Table 1
Configuration of the numerical integration presented in Section5. Δt is nominal stepsize (days), tint is interval time of output (years), and T is total integration period (million years). The parameter combinations of α (or α′) and k2 from (a) to (f) are common to both in Figures 2, 3, 6, and 7.
α or α′ k 2 m ′ / m 0 Δ t t i n t T (a) α = 0.90 0.900 3.0404326 × 1 0 - 6 2 500 10 (b) α = 0.90 0.400 3.0404326 × 1 0 - 6 2 500 10 (c) α ′ = 0.70 0.600 3.0404326 × 1 0 - 6 2 500 10 (d) α ′ = 0.60 0.265 4.7739597 × 1 0 - 5 20 1000 200 (e) α ′ = 0.60 0.260 4.7739597 × 1 0 - 5 20 1000 200 (f) α ′ = 0.80 0.313 9.5479194 × 1 0 - 6 4 500 30In Table1, mass ratios between perturbing body and central body m′/m0 for (a), (b), and (c) are close to the ratio of the Earth-Moon total mass and the solar mass (~1/314 of the mass ratio between Jupiter and the Sun). The m′/m0 value for (f) is nearly 1/100 of the mass ratio of Jupiter and the Sun. The reason for these small m′/m0 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 Δt for these cases by the same reason. As for (d) and (e), the m′/m0 value is close to 1/20 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 m′/m0 as well as larger stepsize Δt. Here, readers should recall the fact that the influence of the mass ratio m′/m0 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 m′ to pm′ (where p is a positive real number); only the timescale of orbital evolution would change from t to t/p. This is why we can arbitrarily change the mass ratio m′/m0 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 bodya′=5.2025217 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) l′ is set to zero at time t=0. By the definition of CR3BP, the perturber’s eccentricity e′ and inclination i′ are both set to zero. As for the perturbed body’s initial orbital elements, semimajor axis a is automatically determined from the value of α or α′. We set the initial values of eccentricity e and argument of pericenter g of perturbed body along the x- and y-axis on the (ecos⁡g,esin⁡g) plane: we selected e of perturbed body from e=0 to emax with the interval of 0.01. As for g, we used four values of g=0, π/2, π, and 3π/2. The values of the perturbed body’s i are automatically determined by relation (53). Initial values of the perturbed body’s longitude of ascending node (h) 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 is(56)d2rdt2+μrr3=∇R,with μ=Gm0. For the outer case we just replace R for R′. 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 Δt 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 codeSWIFT [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 G=0.01720209895, 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 Figure3 whose parameters are described in Table 1. The perturbed body’s (ecos⁡g,esin⁡g) values are plotted on the panels by black and green dots (see below for the reason for this distinction) with the fixed output interval tint designated in Table 1. We chose tint 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 (x,y)=(ecos⁡g,esin⁡g), the ranges of x- and y-axis, meaning of the black dashed circles, and that of the red and blue partial circles.Figure 3
Plots of the perturbed body’s(ecos⁡g,esin⁡g) values obtained from the direct numerical integration whose parameters are listed in Table 1. As for the meaning of the black dashed circles and of the red and blue partial circles, see the caption of Figure 2. The green dots represent the trajectories of bodies that experienced orbital instability during the integration period. Stable trajectories are drawn by black dots.When inspecting Figure3, 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 T. 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 a does not exhibit any major changes throughout the integration period T. Otherwise we categorize the orbital motion as “unstable.” This categorization is done by the naked eye, plotting the time sequence of a for all the perturbed bodies used in the integration (320 bodies per panel, 1920 bodies in total). Browsing through all the plots of a, 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 thex- and y-axis as g=0,π/2,π,3π/2 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 (e,g)=(0.03,0) and (e,g)=(0.03,π) generate almost identical trajectories. For eliminating redundancy in figure and reducing the file size, we did not draw the trajectory from the initial condition (e,g)=(0.03,π) and just drew that from (e,g)=(0.03,0). 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 Figures3 and 2, we say that their agreement is overall good. In Figure 3(a) for the inner CR3BP, existence of stationary points at g=±π/2 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 g=0,π 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 Figure3(b), existence of stationary points at g=±π/2 as well as those at g=0,π is confirmed by the numerical integration. There must be a saddle point at the origin (x,y)=(0,0), 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 Figure3(c) for the outer CR3BP, existence of a stationary point at the origin (x,y)=(0,0) as well as those at g=0,π 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 y-axis.In Figures3(d) and 3(e), the existence of stationary points is again confirmed by the numerical integration. We see a pair of stationary points at g=±π/2 in both panels. The origin (x,y)=(0,0) in Figure 3(d) makes a saddle point, as suggested by the quadrature (Figure 2(d)). The origin (0,0) 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 x-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 (0,0) 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 Figures3(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 Figure3(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 g=±π/2 in Figure 3(f) that we do not see in Figure 2(f). Also, although we see stationary points along the x-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 (e≤emax) 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)).Figure 4
(a) Extended plot of Figure2(f) up to e=emax (numerical quadrature). (b) Extended plot of Figure 3(f) (numerical integration). Note that we set the contour interval in (a) ten times larger than in Figure 2(f) for a better visibility of plots. The interval of initial eccentricity in (b) is twice larger than in Figure 3(f) for the same reason.
(a) (b)So far we have not completely figured out the cause of the difference in trajectory patterns between Figures2(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 5λ′-7λ+2g, 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 5λ′-7λ+2g 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]).Figure 5
Some of the peculiar periodic orbits seen in Figure4(b). The initial e and g are as follows: (a) (e,g)=(0.06,0), (b) (e,g)=(0.08,0), (c) (e,g)=(0.10,0), (d) (e,g)=(0.07,π/2), and (e) (e,g)=(0.36,π/2). The left column: the usual (ecos⁡g,esin⁡g) plots. As for the meaning of the red and blue partial circles, see the caption of Figure 2. The middle column: the argument 5λ′-7λ+2g plotted with a time interval of 1500 years. The right column: osculating semimajor axis a plotted with an interval of 1500 years. All the orbital elements are calculated on the heliocentric coordinate centered on the central mass.
(a) (b) (c) (d) (e)Although all the five perturbed bodies whose trajectories are shown in Figure5 started with the same initial semimajor axis (a0=6.503152AU), we see a systematic difference in the average value of the semimajor axis a between 5(a–c) and 5(d, e) (see the panels in the right column of Figure 5). The short-term oscillation of a 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 a 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 Figure4(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(sin⁡icos⁡h,sin⁡icos⁡h) of the perturbed body where i is its inclination and h 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 cos⁡i takes place when e=0 as cos⁡imin=k  (see Section 4), it is clear that the theoretical largest value of sin⁡i is sin⁡imax=1-k2 when e=0 (note that this value happens to be the same as emax). We drew the outer boundary sin⁡imax=1-k2 by the black dashed circles in Figure 6.Figure 6
Plots of the perturbed body’s(sin⁡icos⁡h,sin⁡isin⁡h) values obtained from the direct numerical integration whose parameters are listed in Table 1. The black dashed circles represent the theoretical maximum sin⁡i value (which is the same as emax). Meanings of the black dots and the green dots remain the same as in Figure 3.In the doubly averaged CR3BP system that we now consider, eccentricitye and inclination i are correlated with each other through the definition of k2 in (53). However, since longitude of ascending node h is not included in the disturbing function R,h does not have a particular correlation with eccentricity e, inclination i, or argument of pericenter g. 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 (sin⁡icos⁡h,sin⁡icos⁡h) plane; i and h are not correlated, and these trajectories do not compose equipotential contours.For more demonstrations of the time variation ofi and h, we picked one body per each of the panels from Figure 6(a) to Figure 6(f) and showed the time variation of their e, i, g, and h in Figure 7. Initial values of (e,g) 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 g=π/2 or 0 (a, b, c, d, and f), but we chose a body whose argument of pericenter circulates from 0 to 2π (e).Figure 7
Variation ofe, i, g, and h of representative bodies in Figures 6(a)–6(f). Initial values of e and g for each body are as follows: (a) (e,g)=(0.14,π/2), (b) (e,g)=(0.50,π/2), (c) (e,g)=(0.34,0), (d) (e,g)=(0.15,π/2), (e) (e,g)=(0.20,0), and (f) (e,g)=(0.06,0). Note that the unit of g and h is degree, not radian, for consistency with the unit of i. All the elements are plotted with a time interval of 1500 years and are calculated on the heliocentric coordinate centered on the central mass.
(a) (b) (c) (d) (e)  (f)The panels in Figure7 immediately tell us the correlation between eccentricity e and inclination i, which is an obvious outcome of the preservation of k2 in (53). The argument of pericenter g 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 h has no correlation with other elements; it shows a monotonic decrease in all the panels.We do not show the analytic counterpart of Figure6 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 h. For example, in Lidov [2], the author uses the mutual dependence between e, i, and g and constructs a definite integral for quadrature that formally depicts the function form of h depending on time. Incidentally, if we limit ourselves to the doubly averaged inner CR3BP at the quadruple-order analytic approximation (2n=2; see Section 3.2), analytic solution of h and that of g (hence those of e and i) using elliptic functions have already been derived in much of the literature (e.g., [3, 5–7, 61, 62]).


## 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 Figure2, 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.Table 1
Configuration of the numerical integration presented in Section5. Δt is nominal stepsize (days), tint is interval time of output (years), and T is total integration period (million years). The parameter combinations of α (or α′) and k2 from (a) to (f) are common to both in Figures 2, 3, 6, and 7.
α or α′ k 2 m ′ / m 0 Δ t t i n t T (a) α = 0.90 0.900 3.0404326 × 1 0 - 6 2 500 10 (b) α = 0.90 0.400 3.0404326 × 1 0 - 6 2 500 10 (c) α ′ = 0.70 0.600 3.0404326 × 1 0 - 6 2 500 10 (d) α ′ = 0.60 0.265 4.7739597 × 1 0 - 5 20 1000 200 (e) α ′ = 0.60 0.260 4.7739597 × 1 0 - 5 20 1000 200 (f) α ′ = 0.80 0.313 9.5479194 × 1 0 - 6 4 500 30In Table1, mass ratios between perturbing body and central body m′/m0 for (a), (b), and (c) are close to the ratio of the Earth-Moon total mass and the solar mass (~1/314 of the mass ratio between Jupiter and the Sun). The m′/m0 value for (f) is nearly 1/100 of the mass ratio of Jupiter and the Sun. The reason for these small m′/m0 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 Δt for these cases by the same reason. As for (d) and (e), the m′/m0 value is close to 1/20 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 m′/m0 as well as larger stepsize Δt. Here, readers should recall the fact that the influence of the mass ratio m′/m0 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 m′ to pm′ (where p is a positive real number); only the timescale of orbital evolution would change from t to t/p. This is why we can arbitrarily change the mass ratio m′/m0 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 bodya′=5.2025217 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) l′ is set to zero at time t=0. By the definition of CR3BP, the perturber’s eccentricity e′ and inclination i′ are both set to zero. As for the perturbed body’s initial orbital elements, semimajor axis a is automatically determined from the value of α or α′. We set the initial values of eccentricity e and argument of pericenter g of perturbed body along the x- and y-axis on the (ecos⁡g,esin⁡g) plane: we selected e of perturbed body from e=0 to emax with the interval of 0.01. As for g, we used four values of g=0, π/2, π, and 3π/2. The values of the perturbed body’s i are automatically determined by relation (53). Initial values of the perturbed body’s longitude of ascending node (h) 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 is(56)d2rdt2+μrr3=∇R,with μ=Gm0. For the outer case we just replace R for R′. 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 Δt 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 codeSWIFT [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 G=0.01720209895, 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 Figure3 whose parameters are described in Table 1. The perturbed body’s (ecos⁡g,esin⁡g) values are plotted on the panels by black and green dots (see below for the reason for this distinction) with the fixed output interval tint designated in Table 1. We chose tint 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 (x,y)=(ecos⁡g,esin⁡g), the ranges of x- and y-axis, meaning of the black dashed circles, and that of the red and blue partial circles.Figure 3
Plots of the perturbed body’s(ecos⁡g,esin⁡g) values obtained from the direct numerical integration whose parameters are listed in Table 1. As for the meaning of the black dashed circles and of the red and blue partial circles, see the caption of Figure 2. The green dots represent the trajectories of bodies that experienced orbital instability during the integration period. Stable trajectories are drawn by black dots.When inspecting Figure3, 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 T. 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 a does not exhibit any major changes throughout the integration period T. Otherwise we categorize the orbital motion as “unstable.” This categorization is done by the naked eye, plotting the time sequence of a for all the perturbed bodies used in the integration (320 bodies per panel, 1920 bodies in total). Browsing through all the plots of a, 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 thex- and y-axis as g=0,π/2,π,3π/2 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 (e,g)=(0.03,0) and (e,g)=(0.03,π) generate almost identical trajectories. For eliminating redundancy in figure and reducing the file size, we did not draw the trajectory from the initial condition (e,g)=(0.03,π) and just drew that from (e,g)=(0.03,0). 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 Figures3 and 2, we say that their agreement is overall good. In Figure 3(a) for the inner CR3BP, existence of stationary points at g=±π/2 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 g=0,π 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 Figure3(b), existence of stationary points at g=±π/2 as well as those at g=0,π is confirmed by the numerical integration. There must be a saddle point at the origin (x,y)=(0,0), 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 Figure3(c) for the outer CR3BP, existence of a stationary point at the origin (x,y)=(0,0) as well as those at g=0,π 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 y-axis.In Figures3(d) and 3(e), the existence of stationary points is again confirmed by the numerical integration. We see a pair of stationary points at g=±π/2 in both panels. The origin (x,y)=(0,0) in Figure 3(d) makes a saddle point, as suggested by the quadrature (Figure 2(d)). The origin (0,0) 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 x-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 (0,0) 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 Figures3(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 Figure3(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 g=±π/2 in Figure 3(f) that we do not see in Figure 2(f). Also, although we see stationary points along the x-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 (e≤emax) 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)).Figure 4
(a) Extended plot of Figure2(f) up to e=emax (numerical quadrature). (b) Extended plot of Figure 3(f) (numerical integration). Note that we set the contour interval in (a) ten times larger than in Figure 2(f) for a better visibility of plots. The interval of initial eccentricity in (b) is twice larger than in Figure 3(f) for the same reason.
(a) (b)So far we have not completely figured out the cause of the difference in trajectory patterns between Figures2(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 5λ′-7λ+2g, 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 5λ′-7λ+2g 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]).Figure 5
Some of the peculiar periodic orbits seen in Figure4(b). The initial e and g are as follows: (a) (e,g)=(0.06,0), (b) (e,g)=(0.08,0), (c) (e,g)=(0.10,0), (d) (e,g)=(0.07,π/2), and (e) (e,g)=(0.36,π/2). The left column: the usual (ecos⁡g,esin⁡g) plots. As for the meaning of the red and blue partial circles, see the caption of Figure 2. The middle column: the argument 5λ′-7λ+2g plotted with a time interval of 1500 years. The right column: osculating semimajor axis a plotted with an interval of 1500 years. All the orbital elements are calculated on the heliocentric coordinate centered on the central mass.
(a) (b) (c) (d) (e)Although all the five perturbed bodies whose trajectories are shown in Figure5 started with the same initial semimajor axis (a0=6.503152AU), we see a systematic difference in the average value of the semimajor axis a between 5(a–c) and 5(d, e) (see the panels in the right column of Figure 5). The short-term oscillation of a 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 a 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 Figure4(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(sin⁡icos⁡h,sin⁡icos⁡h) of the perturbed body where i is its inclination and h 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 cos⁡i takes place when e=0 as cos⁡imin=k  (see Section 4), it is clear that the theoretical largest value of sin⁡i is sin⁡imax=1-k2 when e=0 (note that this value happens to be the same as emax). We drew the outer boundary sin⁡imax=1-k2 by the black dashed circles in Figure 6.Figure 6
Plots of the perturbed body’s(sin⁡icos⁡h,sin⁡isin⁡h) values obtained from the direct numerical integration whose parameters are listed in Table 1. The black dashed circles represent the theoretical maximum sin⁡i value (which is the same as emax). Meanings of the black dots and the green dots remain the same as in Figure 3.In the doubly averaged CR3BP system that we now consider, eccentricitye and inclination i are correlated with each other through the definition of k2 in (53). However, since longitude of ascending node h is not included in the disturbing function R,h does not have a particular correlation with eccentricity e, inclination i, or argument of pericenter g. 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 (sin⁡icos⁡h,sin⁡icos⁡h) plane; i and h are not correlated, and these trajectories do not compose equipotential contours.For more demonstrations of the time variation ofi and h, we picked one body per each of the panels from Figure 6(a) to Figure 6(f) and showed the time variation of their e, i, g, and h in Figure 7. Initial values of (e,g) 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 g=π/2 or 0 (a, b, c, d, and f), but we chose a body whose argument of pericenter circulates from 0 to 2π (e).Figure 7
Variation ofe, i, g, and h of representative bodies in Figures 6(a)–6(f). Initial values of e and g for each body are as follows: (a) (e,g)=(0.14,π/2), (b) (e,g)=(0.50,π/2), (c) (e,g)=(0.34,0), (d) (e,g)=(0.15,π/2), (e) (e,g)=(0.20,0), and (f) (e,g)=(0.06,0). Note that the unit of g and h is degree, not radian, for consistency with the unit of i. All the elements are plotted with a time interval of 1500 years and are calculated on the heliocentric coordinate centered on the central mass.
(a) (b) (c) (d) (e)  (f)The panels in Figure7 immediately tell us the correlation between eccentricity e and inclination i, which is an obvious outcome of the preservation of k2 in (53). The argument of pericenter g 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 h has no correlation with other elements; it shows a monotonic decrease in all the panels.We do not show the analytic counterpart of Figure6 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 h. For example, in Lidov [2], the author uses the mutual dependence between e, i, and g and constructs a definite integral for quadrature that formally depicts the function form of h depending on time. Incidentally, if we limit ourselves to the doubly averaged inner CR3BP at the quadruple-order analytic approximation (2n=2; see Section 3.2), analytic solution of h and that of g (hence those of e and i) 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 Sections4 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 1-e2-2n+1/2 in (24) substantially enhances the magnitude of high-order terms in the analytically expanded disturbing function when eccentricity e is large. For example, when e=0.5 the factor becomes ~8.65 (2n=8), ~27.3 (2n=12), and ~48.6 (2n=14). When e=0.8, the factor reaches ~2127 (2n=8), ~126626 (2n=12), and ~977049 (2n=14), possibly causing artificial oscillation of the higher-order analytic terms in the large-e 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.



---

*Source: 8945090-2016-12-06.xml*