The considered three-body system consists of a massless particle in a close orbit around the central body orbited by a perturbing body in a wider inclined and eccentric orbit (see Figure 1). In this part, the perturbing and the perturbed bodies are assumed to be point masses (the perturbing body in Figure 1 is depicted as an ellipsoid for later purposes), and the central body is an oblate spheroid. We choose the central body’s equatorial plane as the reference plane (the mutual orbit plane of the two massive bodies undergoes secular evolution when their non-spherical gravity is included and is therefore not suitable for the reference plane). The inertial reference frame x , y , z is established with the origin located at the mass center of the central body and with the x − y plane lying in the reference plane, the x -axis initially along the nodal line of the perturbing body’s orbit, and the z-axis along the pole of the central body. For convenient, we represent the perturbing body’s orbital coordinate system as e ^ B , e ^ ⊥ B , h ^ B (the superscript denotes unit vector) where e B is the orbit eccenrticity vector pointing to the periapse, h B is the normalized orbit angular momentum vector, and e ^ ⊥ B complete the righ-hand coordinate system. If we ignore the evolution of the perturbing body’s orbit, this frame is also inertially fixed. Similarly, the particle’s orbital coordinate system can be denoted as e ^ , e ^ ⊥ , h ^ .

The perturbation effects of the perturbing body’s gravity on the particle can be stated in the potential form (called the disturbing function) as R P B = μ B 1 r − r B − r B ⋅ r r B 3 (1) where μ B = G m B is the perturbing body’s gravitational parameter; r B and r are the position vectors of the perturbing body and the particle with respect to the central body’s mass center, respectively.

In hierarchical systems, the disturbing function can be expanded into an infinite series in the form of Legendre polynomials as R P B = μ B r B ∑ n = 2 ∞ r r B n P n r ^ ⋅ r ^ B (2) where P n x are the Legendre polynomials about x . According to this formula, the disturbing function can be approximated by truncating the expansion to an arbitrary order, among which the lowest-order approximation, i.e., the quadrupole order, is most widely used in previous studies. However, if the three-body system is not highly hierarchical, especially when considering an eccentric perturbing body, the next several higher orders may be surprisingly important. To this end, here, we truncate the expansion to the hexadecapole order as R P B = − μ B 2 r B { r r B 2 3 r ^ B ⋅ r ^ 2 − 1 + r r B 3 5 r ^ B ⋅ r ^ 3 − 3 r ^ B ⋅ r ^ + r r B 4 35 r ^ B ⋅ r ^ 4 − 30 r ^ B ⋅ r ^ 2 + 3 } (3)

The disturbing function of the oblateness of the central body can be represented by the gravitational coefficient J 2 as R J 2 = μ A J 2 r E 2 2 r 3 1 − 3 r ^ ⋅ z 2 (4) where μ A = G m A is the gravitational parameter of the central body and r E is its mean radius.

Then, the full Newtonian equations of motion, including perturbations of the central body’s oblateness and the third-body gravity, written in the inertial reference frame, are given as r ¨ = − μ A r 3 r + ∂ R J 2 ∂ r + ∂ R P B ∂ r (5) where the disturbing functions R P B and R J 2 are given in Eqs. 1 and 4. If we replace the R P B in Eq. 5 with that in Eq. 3, we can obtain the hexadecapole-order (non-averaged) equations.

In order to reduce the degrees of the system and provide analytical insight into the perturbation effects, the averaging method is proposed to eliminate the short-period terms and derive the equations that capture the long-term evolution of the system. According to the Krylov–Bogoliubov–Mitropolski method of averaging (Bogoliubov and Mitropolsky 1961; Rosengren, 2014), the averaging of the disturbing function is defined as R ¯ α = 1 2 π ∫ 0 2 π R α , M d M (6) where α denotes the orbital elements except the mean anomaly M , which is proportional to time.

The validity of the averaging requires that α is nearly constant in the integration process. Physically, it requires the perturbation forces to be sufficiently small so that it follows a Keplerian orbit within one orbital period. In addition, in the current system, there are two frequencies related to the inner and outer orbital motion respectively. Then, it also requires that the two frequencies are away from each other and there is no resonance between them (the averaged model is not applicable when mean motion resonance occurs). The mean motions of the inner and outer orbits are given as n = μ A / a 3 (7) N B = μ A B / a B 3 (8) where a is the inner orbit semimajor axis; μ A B = G m A + m B ; a B is the perturbing body’s semimajor axis. In the previous studies treating highly hierarchical systems, the conditions of the averaging are rarely discussed in detail. However, considering a stronger perturbed system with moderate hierarchy, these will be no longer negligible. In dealing with secular dynamics of a orbiter of an binary asteroid system, Wang and Fu (2020) have investigated the specific conditions of the averaging for systems with combined effects of the central body’s oblateness and third-body gravity. They concluded that, approximately, the perturbing forces should not exceed 1 per cent of the dominated central body’s gravity, and n / N B should be less than 0.4. Therefore, for a given three-body system possessing the geometry in Figure 1, one can easily determine the range of the orbital altitude that ensures the reliability of the averaging. In the following calculations, we assume the averaging is performed within the valid range of the orbital altitude.

By the averaging theory, the secular Lagrange planetary equations in terms of vector elements, taking the compact and symmetrical form, are given by h ˙ = h × ∂ R ¯ ′ ∂ h + e × ∂ R ¯ ′ ∂ e e ˙ = e × ∂ R ¯ ′ ∂ h + h × ∂ R ¯ ′ ∂ e (9) where R ¯ ′ is the averaged disturbing function scaled by μ A a (where a is constant in the averaged dynamics); h = H / μ A a is the scaled orbit angular momentum; e is the eccentricity vector. h and e in Eq. 9 denote mean elements (we drop the over bar for simplicity).

In the hierarchical three-body systems, as stated above, there are two timescales associated with the orbital motion n and N B respectively. Under certain conditions, averaging over the two timescales is feasible. We perform the first averaging over the mean anomaly of the inner orbit. The quadrupole order in Eq. 3 is averaged as R ¯ Q u a = − μ B 2 r B 3 3 r ^ B ⋅ r r ¯ ⋅ r ^ B − r 2 ¯ (10) where the product of arbitrary two vectors a b is called a dyad and, in the language of the matrix algebra, is equivalent to a b T . The related terms are averaged as (Rosengren and Scheeres, 2013) r 2 ¯ = a 2 1 + 3 2 e 2 (11) r r ¯ = 1 2 a 2 1 + 4 e 2 e ^ e ^ + 1 − e 2 e ^ ⊥ e ^ ⊥ (12)

Substituting Eqs. 11–12 into Eq. 10, disregarding the irrelevant constant terms and scaling to μ A a , we have R ¯ Q u a ′ = 3 μ B 4 n r B 3 5 r ^ B ⋅ e 2 − r ^ B ⋅ h 2 − 2 e 2 (13) by virtue of the dyadic identity e ^ e ^ + e ^ ⊥ e ^ ⊥ = I − h ^ h ^ in which I denotes the identity matrix.

Averaging of the octupole-order terms in Eq. 3 requires (Rosengren, 2014) r ¯ = − 3 2 a e (14) r 2 r ¯ = − 5 8 a 3 4 + 3 e 2 e (15) r r r ¯ = − 5 8 a 3 e 3 + 4 e 2 e ^ e ^ e ^ + 1 − e 2 e ^ ⊥ e ^ ⊥ e ^ + e ^ e ^ ⊥ e ^ ⊥ + e ^ ⊥ e ^ e ^ ⊥ (16) where the product of three vectors a b c is called a triad. To illustrate how it is calculated with other vectors, considering the averaging of r ^ B ⋅ r 3 and by Eq. 16, we have r ^ B ⋅ r 3 ¯ = − 5 8 a 3 7 r ^ B ⋅ e 3 + 3 r ^ B ⋅ e h 2 − r ^ B ⋅ h 2 (17) by dotting r ^ B with each of the vector elements in each triad on the right-hand side of Eq. 16. Through Eqs. 14–17, the octupole order is averaged as R ¯ O c t ′ = − 5 μ B a 16 n r B 4 35 r ^ B ⋅ e 3 + 3 1 − 8 e 2 r ^ B ⋅ e − 15 r ^ B ⋅ e r ^ B ⋅ h 2 (18)

The averaging of the hexadecapole order follows a similar procedure and is elaborated in Appendix B. Then, combining Eqs. 10 and 18, and (B5), the single-averaged third-body disturbing function to the hexadecapole order is given as R ¯ P B ′ = R ¯ Q u a ′ + R ¯ O c t ′ + R ¯ H e x ′ (19)

Performing the averaging over the disturbing function of the central body’s oblateness, Eq. 4, gives R ¯ J 2 = μ A J 2 r E 2 2 3 z ⋅ r ^ r ^ r 3 ¯ ⋅ z − 1 r 3 ¯ (20) in which, we note 1 r 3 ¯ = 1 a 3 h 3 (21) r ^ r ^ r 3 ¯ = 1 2 a 3 h 3 I − h ^ h ^ (22) Then, yields R ¯ J 2 ′ = n J 2 r E 2 4 a 2 h 3 3 z ⋅ h ^ 2 − 1 (23) Finally, by substituting Eqs. 19–23 into the secular Lagrange planetary equations, Eq. 9, the single-averaged equations of motion in terms of h and e , with the combined effects of the third-body gravity and the central body’s oblateness, can be obtained (we do not show the explicit expressions here and please refer to Wang and Fu (2020) if necessary). The single-averaged dynamical model is derived by only averaging over the orbital motion of the inner orbit. That is, this model provides chance to study the resonance between the secular precession of the inner orbit and the mean motion of the outer body. On the other hand, the single-averaged model may be significantly more accurate than the double-averaged model in certain conditions (Wang and Gurfil, 2016; Nie et al., 2019; Carvalho et al., 2022).

Here, we do not discuss the applications of the single-averaged model and directly perform the second averaging over the outer orbit. Since the inclination and eccentricity of the perturbing body are included, the second-averaged disturbing functions would be different from previous studies that assuming a coplanar and circular configuration (Rosengren and Scheeres, 2013; Fu and Wang, 2021). The eccentricity e B and inclination i B are implied in the vectorial elements e B and h B .

The second averaging of the quadrupole-order approximation, Eq. 13, is given as R = Q u a ′ = 3 μ B 4 n 5 e ⋅ r ^ B r ^ B r B 3 ¯ ⋅ e − h ⋅ r ^ B r ^ B r B 3 ¯ ⋅ h − 2 e 2 1 r B 3 ¯ (24) in which we note 1 r B 3 ¯ = 1 a B 3 h B 3 (25) r ^ B r ^ B r B 3 ¯ = 1 2 a B 3 h B 3 I − h ^ B h ^ B (26) and thus, we have R = Q u a ′ = 3 μ B 8 n 1 a B 3 h B 3 e 2 − 5 e ⋅ h ^ B 2 − h 2 + h ⋅ h ^ B 2 (27)

This is the classical double-averaged third-body disturbing function in the Hill approximation in vectorial elements. Substituting it into the Lagrange planetary equations, Eq. 9, the quadrupole-level equations of motion are derived as h ˙ Q u a = 3 μ B 4 n 1 a B 3 h B 3 5 e ⋅ h ^ B h ^ B × e − h ⋅ h ^ B h ^ B × h (28) e ˙ Q u a = 3 μ B 4 n 1 a B 3 h B 3 2 h × e + 5 e ⋅ h ^ B h ^ B × h − h ⋅ h ^ B h ^ B × e (29) Note that, although we have included the eccentricity of perturbing body’s orbit, there is no components along e ^ B or e ^ ⊥ B , which yields a dynamical axisymmetry about the pole of the perturbing body’s orbital plane (i.e., h ^ B ), leading to the conservation of the z -component of the inner orbit’s angular momentum.

The second averaging over the octupole order, Eq. 18, requires r ^ B r B 4 ¯ = e B a B 4 h B 5 (30) r ^ B r ^ B r ^ B r B 4 ¯ = e B 4 a B 4 h B 5 3 e ^ B e ^ B e ^ B + e ^ ⊥ B e ^ ⊥ B e ^ B + e ^ B e ^ ⊥ B e ^ ⊥ B + e ^ ⊥ B e ^ B e ^ ⊥ B (31) and through which, we have R = O c t ′ = − 5 μ B a 64 n e B a B 4 h B 5 105 e ⋅ e ^ B e 2 − e ⋅ h ^ B 2 + 12 1 − 8 e 2 e ⋅ e ^ B + 30 e ⋅ h ^ B h ⋅ e ^ B h ⋅ h ^ B − 15 e ⋅ e ^ B h 2 − h ⋅ h ^ B 2 (32) The secular equations of motion in the octupole order are thus given as h ˙ O c t = − 15 μ B a 64 n e B a B 4 h B 5 − 10 e ⋅ h ^ B h ⋅ h ^ B e ^ B × h − 10 e ⋅ h ^ B h ⋅ e ^ B + e ⋅ e ^ B h ⋅ h ^ B h ^ B × h + 35 e ⋅ h ^ B 2 − 5 h ⋅ h ^ B 2 − 8 e 2 + 1 e ^ B × e + 10 7 e ⋅ e ^ B e ⋅ h ^ B − h ⋅ e ^ B h ⋅ h ^ B h ^ B × e (33) e ˙ O c t = − 15 μ B a 64 n e B a B 4 h B 5 − 10 e ⋅ h ^ B h ⋅ h ^ B e ^ B × e − 10 e ⋅ h ^ B h ⋅ e ^ B + e ⋅ e ^ B h ⋅ h ^ B h ^ B × e + 35 e ⋅ h ^ B 2 − 5 h ⋅ h ^ B 2 − 8 e 2 + 1 e ^ B × h + 10 7 e ⋅ e ^ B e ⋅ h ^ B − h ⋅ e ^ B h ⋅ h ^ B h ^ B × h + 16 e ⋅ e ^ B h × e (34) Note that if e B = 0 , the octupole order vanishes. When e B > 0 , the components along e ^ B violate the axisymmetry of the third-body perturbation predicted by the quadrupole-level model. As a result, the system is no longer integrable and the octupole order qualitatively changes the secular evolutionary behaviors of the system.

The second averaging and the corresponding secular equations of the hexadecapole order are presented in Supplementary Appendix B. As we can see, they also contain components in the e ^ B direction, but when e B = 0 , these terms disappear and the axisymmetry is retrieved. It turns out that the axisymmetry only holds in the limit of a massless particle perturbed by a third body on a circular orbit (Naoz et al., 2013).

Finally, applying the averaged disturbing function of the oblateness to the secular equations, we obtain h ˙ J 2 = − 3 n J 2 r E 2 2 a 2 h 5 z ⋅ h z × h (35) e ˙ J 2 = 3 n J 2 r E 2 4 a 2 h 5 1 − 5 h 2 z ⋅ h 2 h × e + 2 z ⋅ h z × e (36) So far, we have obtained the complete double-averaged equations of motion, including the secular effects of the central body’s oblateness and the third-body gravity (to the hexadecapole order), in the form of vector elements, as h ˙ = h ˙ J 2 + h ˙ Q u a + h ˙ O c t + h ˙ H e x e ˙ = e ˙ J 2 + e ˙ Q u a + e ˙ O c t + h ˙ H e x (37) In order to simplify the analysis of the perturbing effects of different perturbing forces, we rewrite the resultant double-averaged disturbing function as R = ′ = ε J 2 R ¯ J 2 * + ε P B R = Q u a * + ε O c t R = O c t * + ε H e x R = H e x * (38) where ε J 2 and ε P B are the coefficients of the oblateness and the quadrupole-level third-body gravity respectively, defined as ε J 2 = n J 2 r E 2 a 2 (39) ε P B = μ B n 1 a B 3 h B 3 (40) ε O c t and ε H e x are the octupole and hexadecapole order coefficients respectively, defined as ratio of the octupole and hexadecapole order to the quadrupole order ε O c t = ε O c t ε P B = a a B e B h B 2 (41) ε H e x = ε H e x ε P B = a 2 a B 2 1 1 − e B 2 2 (42) To eliminate an extra parameter, we divide the disturbing function by the coefficient ε P B , R = = κ R ¯ J 2 * + R = Q u a * + ε O c t R = O c t * + ε H e x R = H e x * (43) where κ is defined as the ratio of the oblateness perturbation to the third-body perturbation κ = ε J 2 ε P B = m A m B J 2 r E 2 a 2 a B 3 h B 3 a 3 (44) The equations of motion for the disturbing function Eq. 43 is equivalent to apply a new rescaled time to Eq. 37. The rescaled time is defined as τ = t × ε P B (45) where t is the true time. As a consequence, the double-averaged equations of motion with respect to the rescaled time τ are written as h ˙ = κ h ˙ J 2 * + h ˙ Q u a * + ε O c t h ˙ O c t * + ε H e x h ˙ H e x * e ˙ = κ e ˙ J 2 * + e ˙ Q u a * + ε O c t e ˙ O c t * + ε H e x e ˙ H e x * (46) The coefficients κ and ε O c t and the parameter i B can be used to characterize the secular evolution of the system (the hexadecapole order ε H e x typically has smaller effects). (1) For the special cases with i B = 0 , if we take κ = 0 ( ε J 2 = 0 , namely, only the third-body perturbation included), it degenerates into the classical Lidov-kozai mechanism ( ε O c t ∼ 0 ) or the eccentric Lidov-Kozai effects ( ε O c t > 0 ); if we take κ → ∞ ( ε P B = 0 , i.e., only the oblateness effects considered), we have the well-known J 2 -perturbed problem. Delsate et al. (2010) has conducted a global analysis on the phase structures of the secular dynamics of the integrable case where J 2 and the third-body perturbation are included but assuming i B ∼ ε O c t ∼ 0 . On the other hand, it indicates that, when ε O c t > 0.01 , the octupole order is non-negligible and typically leads to qualitatively new features such as orbit flips (Lithwick and Naoz, 2011). (2) When i B > 0 , the misalignment between the central body’s pole and the perturbing body’s orbit angular momentum also induces new behaviors, but in a much different way (Yokoyama, 1999). That is, both ε O c t and i B qualitatively change the global behavior of the mentioned integrable case, and it would be interesting to view their combined secular effects.

Observations indicate that small binary asteroid systems typically contain a fast-spinning, flattened primary and a synchronously rotating, elongated secondary (Ćuk and Nesvorný, 2010). This is the natural outcome of the rotation-fission formation theory based on the Yarkovsky-Radzievskii-Paddack effect (YORP) and the evolution mechanism due to the mutual body tides (Jacobson, 2011). Therefore, it is feasible to make some assumptions about these systems so as to provide analytical insights into the secular dynamics: (1) modeling the primary as an oblate spheroid and the secondary as a triaxial ellipsoid, both rotating around their maximum inertial axis; (2) assuming the secondary is in a 1:1 spin-orbit resonance and the primary lying along its minimum inertial axis, i.e., the long-axis configuration (Scheeres, 2013); and (3) considering the secondary in an inclined and eccentric orbit around the primary (Figure 1). This similar model has been employed in many previous studies in investigating the general dynamics of a binary asteroid system as well as an artificial orbiter, with the massive two bodies restricted to planar motion (Bellerose and Scheeres, 2008; Chappaz and Howell, 2015; Dell’Elce et al., 2017; Wang and Hou, 2020; Wang and Fu, 2020; Fu and Wang, 2021). However, the non-planar oblate spheroid-ellipsoid configuration, i.e., including the inclination of the mutual orbit relative to the primary’s equator, has seldom been focused.

One remark is that when the synchronously rotating secondary is in an eccentric orbit, there is a critical orbit eccentricity and a libration amplitude of the secondary beyond which the 1:1 spin-orbit resonance is unstable (Wang and Hou, 2020). In the present paper, we consider the mutual orbit eccentricity lower than the critical value but ignore the libration angle. This is feasible because, as we will show later, the non-spherical gravity terms of the secondary appear as “hexadecapole order” (Eq. 57), and the effect of the libration angle is believed to be even higher order and thus can be neglected.

On the other hand, although evidence implies that synchronous secondary generally has a small orbit eccentricity (Pravec et al., 2016; Wang and Hou, 2020), the lower hierarchy of these systems causes the octupole order to be significant. That is, the effects of the secondary’s eccentricity are amplified in the binary asteroid systems compared to those in highly hierarchical systems (such as the Orbiter-Mercury-Sun system).

Considering the two massive bodies in an arbitrary configuration (Figure 2), with the axes of their body-fixed frames along their principle axes of inertial, the second degree and order mutual potential is written as (Wang and Fu (2020)) U r B , C = − G m A m B r B + G m A m B 2 r B 3 τ 0 A − 3 τ 0 A r ^ B ⋅ z 2 + τ 0 B − 3 τ 0 B r ^ B ⋅ γ 2 − 6 τ 2 B r ^ B ⋅ α 2 − r ^ B ⋅ β 2 (47) where τ 0 A = C 20 , A R A 2 ( J 2 = − C 20 ), τ 0 B = C 20 , B R B 2 , and τ 2 B = C 22 , B R B 2 , in which R A and R B are the mean equatorial radius, and C 20 and C 22 are the coefficients of the second degree and order gravity; C = α , β , γ is the secondary’s attitude matrix relative to the primary, or, since the rotation of the primary is trivial, we define C as the transform matrix from the secondary’s body-fixed reference frame u , v , w to the primary centered inertial reference frame x , y , z (see Supplementary Appendix A).

Then, in the synchronous configuration where r B and u are aligned (see Figure 1), the mutual potential is reduced to U r B = − G m A m B r B + G m A m B 2 r B 3 τ 0 A − 3 τ 0 A r ^ B ⋅ z 2 + τ 0 B − 6 τ 2 B (48) Since the relative attitude of the two bodies is fixed (we assume that the slow precession of the mutual orbit does not break this equilibrium configuration), the relative motion of the two bodies, under their mutual gravitation, can be modeled as the perturbed Keplerian two-body problem, in which the perturbation is their non-spherical mutual gravity. The relative motion is governed by r ¨ B = − μ A B r B 2 r ^ B + 3 μ A B 2 r B 4 τ 0 A 1 − 5 r ^ B ⋅ z 2 r ^ B + 2 r ^ B ⋅ z z + τ 0 B − 6 τ 2 B r ^ B (49) where μ A B = G m A + m B . Then, we perform averaging over the mutual orbit, and by the Lagrange planetary equations in terms of classical elements, the secular equations of motion are given as a ˙ B = i ˙ B = e ˙ B = 0 (50) Ω ˙ B = 3 n B τ 0 A 2 a B 2 1 − e B 2 2 cos ⁡ i B (51) ω ˙ B = 3 n B 4 a B 2 1 − e B 2 2 τ 0 A 1 − 5 cos ⁡ i B − 2 τ 0 B − 6 τ 2 B τ 0 A (52) As we can see, the shape of the secondary’s orbit is unchanged but its orientation undergoes precession due to the secular evolution of the argument of perigee ω B and longitude of the ascending node Ω B . Note that the precession of Ω B is caused by the primary’s oblateness, while the evolution of ω B is a joint effect of the primary’s and secondary’s non-spherical gravity. Given the initial orbital conditions, Ω B and ω B vary linearly with time, and we can easily add them into the particle’s orbital motion equations.

We consider the orbital dynamics around the primary (an oblate spheroid) with the secondary (an ellipsoid) as the perturbing body. As is known, the third-body perturbing acceleration is the difference between the absolute acceleration that the particle feels and the absolute acceleration that the central body feels a P B = a B − a B − A (53) where a B is given as a B = − μ B ρ 2 ρ ^ + 3 μ B τ 0 B 2 ρ 4 1 − 5 ρ ^ ⋅ w 2 ρ ^ + 2 ρ ^ ⋅ w w − 3 μ B τ 2 B ρ 4 5 ρ ^ ⋅ u 2 − ρ ^ ⋅ v 2 ρ ^ − 2 ρ ^ ⋅ u u − ρ ^ ⋅ v v (54) in which ρ = r − r B ; from Eq. 48, a B − A is given as a B − A = μ B r B 2 r ^ B − 3 μ B 2 r B 4 τ 0 A 1 − 5 r ^ B ⋅ z 2 r ^ B + 2 r ^ B ⋅ z z − 3 μ B 2 r B 4 τ 0 B − 6 τ 2 B r ^ B (55) For use in perturbation analysis, it is convenient to recast a P B as a disturbing function R P B = μ B r − r B − μ B τ 0 B 2 r − r B 3 1 − 3 r ⋅ γ 2 r − r B 2 + 3 μ B τ 2 B r − r B 5 r − r B ⋅ α 2 − r ⋅ β 2 − μ B r B 2 r ^ B ⋅ r + 3 μ B τ 0 A 2 r B 4 1 − 5 r ^ B ⋅ z 2 r ^ B ⋅ r + 2 r ^ B ⋅ z z ⋅ r + 3 μ B 2 r B 4 τ 0 B − 6 τ 2 B r ^ B ⋅ r (56) Substituting R P B in Eq. 56 into Eq. 5, we obtain the extended full equations of motion including the non-spherical third-body gravity.

Then, by using Taylor expansion in terms of r / r B and truncated to the hexadecapole order, the disturbing function is approximated as R P B = − μ B 8 r B 4 r / r B 2 3 r ^ B ⋅ r ^ 2 − 1 + 4 r / r B 3 5 r ^ B ⋅ r ^ 3 − 3 r ^ B ⋅ r ^ + r / r B 4 35 r ^ B ⋅ r ^ 4 − 30 r ^ B ⋅ r ^ 2 + 3 − 3 μ B 2 r B τ 0 A r r B 3 1 − 5 r ^ B ⋅ z 2 r ^ B ⋅ r ^ + 2 r ^ B ⋅ z r ^ ⋅ z + 3 μ B 4 r B τ 0 B r 2 r B 4 5 r ^ B ⋅ r ^ 2 − 2 r ^ ⋅ γ 2 − 1 − 3 μ B 2 r B τ 2 B r 2 r B 4 17 r ^ B ⋅ r ^ 2 − 2 r ^ ⋅ β 2 − 5 (57) The terms in the first curly brace correspond to the point mass gravity, which is also given in Eq. 3. The other terms are the contributions of the non-spherical gravity of the two massive bodies. As is clear, if the two massive bodies are fairly oblate or elongated, the strength of the primary’s oblateness, represented by τ 0 A , may be comparable to the point mass’s octupole order, while the secondary’s second degree and order gravity, τ 0 B and τ 2 B , is comparable to the hexadecapole order. However, in general cases, these non-spherical terms are smaller than other terms of the same corresponding order.

Calculating the averaging of these non-spherical terms is similar to those in Section 2.2, and for brevity, we illustrate them in Supplementary Appendix C. Finally, we obtain the extended double-averaged equations of motion as h ˙ = κ h ˙ J 2 * + h ˙ Q u a * + ε O c t h ˙ O c t * + ε O c t N h ˙ O c t N + ε H e x h ˙ H e x * + ε H e x , 0 N h ˙ H e x , 0 N ∗ + ε H e x , 2 N h ˙ H e x , 2 N ∗ e ˙ = κ e ˙ J 2 * + e ˙ Q u a * + ε O c t e ˙ O c t * + ε O c t N e ˙ O c t N ∗ + ε H e x e ˙ H e x * + ε H e x , 0 N e ˙ H e x , 0 N ∗ + ε H e x , 2 N e ˙ H e x , 2 N ∗ (58) where ε O c t N = − τ 0 A / a 2 × ε O c t , ε H e x , 0 N = − τ 0 B / a 2 × ε H e x , and ε H e x , 2 N = τ 2 B / a 2 × ε H e x are the coefficients of the non-spherical effects of the third-body perturbation.

The applications to the binary asteroid system 4951 Iwamoto are illustrated as an example to validate our double-averaged model and to present numerical results. Observations indicate that 4951 Iwamoto is a doubly synchronous binary asteroid system with a period of 118 h (Pravec et al., 2016). The diameters of the primary and secondary are inferred to be 5.52 km and 3.34 km, respectively, and the secondary is about 31 km from the primary. Since further information is still unknown, some assumptions are made on the physical and orbital characteristics of the binary system, as shown in Table 1. First, we assume the primary is shaped like an oblate spheroid and the secondary is nearly a tri-axial ellipsoid. Second, the two massive bodies have the equilibrium long-axis configuration, i.e., the secondary is synchronously rotating with the primary lies on its long axis. The secondary is assumed to be in an inclined and eccentric orbit around the primary with e B = 0.2 and i B = 10 ° .

The second degree and order gravity coefficients in Table 1 are calculated by (Wang and Fu, 2020) C 20 = 1 5 R e 2 s 2 − p 2 + q 2 2 (59) C 22 = 1 20 R e 2 p 2 − q 2 (60) where s, q, and p are the semi-major axes of an ellipsoid, satisfying s ≤ q ≤ p .

As stated above, strong perturbations may invalidate the averaging process. It means, for the combined effects of the oblateness of the central body and the gravity of the third body, there are lower and upper limits on the orbit altitude that ensure reliability of the averaging. According to Wang and Fu (2020), for the system parameters given in Table 1, the valid range of the orbit semimajor axis around the primary is about 5.0 km–8.0 km, corresponding to the coefficient κ varying from 3.57 to 0.34, as shown in Figure 3. Note that the magnitudes of the octupole order coefficient ε O c t and the hexadecapole order coefficient ε H e x are about the same order (varying between 0.02 and 0.08), due to the octupole order coefficient having a factor of the perturbing body’s eccentricity e B . Therefore, it is necessary to take both of them into the secular model for the purpose of propagating long-term orbital motion.

Comparisons are made among the numerical integrations of the four dynamical models: the full, the quadrupole-level, the octupole-level, and the hexadecapole-level models. The full model (non-averaged) given by Eq. 5 is considered to be exact. The other three models are all double averaged but retained to different levels of approximation.

In Figures 4, 5, the orbital evolution over scaled time is depicted, with the semimajor axes set as 5.3 km and 7.8 km respectively (both within the range where the averaging is reliable). The other initial orbital elements are explicitly given in the caption of Figure 4. Actually, they can be chosen arbitrarily for the purpose of demonstrations of the different models. For comparisons of the performance of these models at different orbital altitudes, the initial conditions in Figure 4; Figure 5,; Figure 6 are kept the same except for the semimajor axis. Note that the true time t is related to the scaled time τ via t = τ / ε P B (Eq. 45). The scaled time of 10 in Figures 4, 5 corresponds to the true time of 330.3 days and 188.6 days, respectively. The precession of the secondary’s orbit is incorporated in the numerical integrations of our full and double-averaged models by including the evolution of Ω B and ω B described in Eqs. 51, 52 respectively. Since the semimajor axis is constant in the averaged dynamics, only four remaining orbital elements: eccentricity e , inclination i , argument of perigee ω , and longitude of ascending node Ω , need to be considered.

From Figures 4, 5, we find that within the valid range of orbital altitude, the hexadecapole-level model could well track the secular orbital evolution, while the quadrupole-level and octupole-level models may have significant deviations. However, beyond this range, for example, as shown in Figure 6, when a = 9.0  km, orbits are strongly perturbed by the third-body gravity and, as a consequence, the large-amplitude short-period oscillation invalidates the averaging processes.

Remember that in the classical Lidov-Kozai mechanism or the extended “ J 2 included Lidov-Kozai mechanism”, the conservation of vertical-component orbit angular momentum, h z , leads to regular, anti-correlated oscillation of the eccentricity and inclination. In the current model, there is a similar anti-correlation behavior between e and i , but not regular. In particular, there is a secular effect, as illustrated in Figure 4A, that may force h z to undergo (chaotic) secular changes, resulting in high eccentricities. Recall that in the eccentric Lidov-Kozai effects, the octupole order can modulate the Lidov-Kozai cycles (the secular change of h z can be considered as the Lidov-Kozai cycles undergo secular drift), and thus excite the extremely high eccentricities and even orbit flips (Katz et al., 2011; Lithwick and Naoz, 2011). However, as shown in Figure 7, our results may be a different secular dynamical effect because they are a joint effort of the third-body eccentricity and inclination. It seems that the large eccentricity increases should be attributed to i B , while e B changes the phase of this secular variation.

In Figure 8, we depict the long-term orbital evolution with different inclinations at a = 5.5  km. For simplicity, only the eccentricity evolution is displayed. First, we show that the hexadecapole-level model is highly accurate, while the quadrupole-level and octupole-level models have moderate deviations. As we have mentioned above, the system is integrable in the coplanar configuration in the quadrupole-level model, even if the perturbing body is eccentric. For comparisons, we integrate these integrable equations which are obtained by substituting i B = 0 into the quadrupole-level model with other parameters remaining unchanged. From Figure 8, we can see the significant differences between the full dynamics and the integrable approximation. In addition, in these integrable cases, the equilibrium points (i.e., frozen orbits) and their stabilities are well-defined. According to Delsate et al. (2010), by using the current system parameters, the analytical unstable inclination interval at a = 5.5  km is calculated as 55.3 ° < i < 81.9 ° and 100.0 ° < i < 124.4 ° , which is symmetric about i = 90 ° . However, in the inclined third-body problem, the inclination is different with respect to the central body’s equator and with respect to the third-body orbit plane, especially when the perturbing body’s orbital plane is in precession, as shown in Figure 8. Thus, if we use the solutions of the integrable cases as initial guesses, the unstable inclination regions should be revised according to the specific geometric configuration of the system.

Then, we show comparisons among these different levels of approximation models for the special cases in which the perturbing body has zero inclination or eccentricity. The corresponding double-averaged equations of motion can be obtained by substituting e B = 0 or i B = 0 into Eq. 58.

Figure 9 illustrates the case when i B = 0 and e B = 0.2 . For the given initial conditions, the orbit is initially located near a stable Lidov-Kozai equilibrium ( ω = 90 ° ) solved by the quadrupole-level model. However, in the full and octupole-level (or higher-level) models, the amplitude of the Lidov-Kozai oscillation gradually increases and finally causes an impact on the primary. In addition, although the octupole order and hexadecapole order have similar magnitudes ε O c t ∼ ε H e x = 0.04 , the dramatic, qualitative changes in dynamics only happen when moving from quadrupole-level to octupole-level models. This is because the octupole order recovers the geometric non-axisymmetry of the system, thus eliminating h z as a constant, leaving only two integrals of motion: the semi-major axis and the double-averaged perturbing potential. That is, there are no other integrals of motion to be eliminated in the higher-level models, and thus no more qualitative changes are expected (Naoz et al., 2013).

By using the same initial conditions, we set e B = 0 and i B = 10 ° in Figure 10. In contrast to the eccentric third-body perturbation problem in Figure 9, when i B ≠ 0 , there are no special integrals of motion in the lowest-order (quadrupole-level) model compared to the higher-level models. Thus, the high order terms may contribute no qualitative changes but quantitatively affect the orbit evolutionary behavior, as shown in Figure 10.

The derived semi-analytical model can also be applied to astrophysical systems covering a wide range of physical scales, such as the exoplanet around an oblate host star with a distant stellar mass or giant planet perturber, where the inner planet can be regarded as massless. Similar dynamical settings have been extensively considered in the eccentric Lidov-Kozai mechanism (Katz et al., 2011; Lithwick and Naoz, 2011; Naoz et al., 2013; Naoz, 2016).

In the scaled dynamical model, Eq. 46, the orbital dynamics can be characterized by three parameters: the oblateness coefficient κ , the octupole-order coefficient ε O c t , and the inclination of the third-body orbit i B . Here, we consider two cases: (1) κ = 0.1 ; ε O c t = 0.01 ,; i B = 10 ° in Figure 11, where ε O c t is much smaller than κ ; and (2) κ = 0.1 , ε O c t = 0.1 , and i B = 10 ° in Figure 12 where ε O c t is comparable to κ .

From Figures 11, 12, we show the important role of the octupole-order terms, especially when the octupole-order coefficient ε O c t is comparable to the oblateness-related coefficient κ (Figure 12), which may lead to the dramatic, qualitatively different evolutionary behaviors. We note that the secular effect shown in Figure 12 is different from that in Figure 7 because the orbit inclination oscillates in a small range even when the orbit reaches high eccentricities in Figure 7, while in Figure 12, the large increases in eccentricities are accompanied by wide excursions of inclination. We suppose that these two secular effects may be the mechanism of forming close-in exoplanets of different inclination types. However, the detailed analysis is beyond the scope of our current work, which will be focused on in our future work.