# Semi-analytical orbital model around an oblate body with an inclined eccentric perturber

- School of Astronautics, Beihang University, Beijing, China

The semi-analytical model (based on the averaging technique) for long-term orbital evolution has proven to be useful in both astrophysical and astrodynamical contexts. In this secular approximation, orbits exhibit rich evolutionary behaviors under the effects of various perturbations. For example, in the hierarchical three-body systems, the Lidov-Kozai mechanism based on the quadrupole-level third-body perturbation model can produce large-amplitude oscillations of the eccentricity and inclination. In recent years, the octupole order has been found to induce dramatically new features when the perturbing body’s orbit is eccentric, including extremely high eccentricities and orbit flips between prograde and retrograde. Motivated by the striking effects of the octupole-order terms, we intend to derive a more general dynamical model by incorporating *J*_{2} of the central body and the inclined eccentric third-body perturbation to the hexadecapole order with its non-spherical gravity also included. This issue can be relevant for astrophysical and astrodynamical systems such as planets in stellar binaries, irregular satellites in planetary systems, and artificial probes about binary asteroid systems. Applications to the binary asteroid system 4951 Iwamoto and a fictitious exoplanetary system are illustrated as examples to validate our secular model. From these numerical results, we show the high accuracy of our secular model. Also, we show the important role of these high-order terms and the effects of the third-body’s inclination and eccentricity. Besides, we have found several different secular effects that could drive large eccentricities.

## 1 Introduction

The secular perturbation theory was initially motivated by studying planetary and planetary satellite motions in the Solar system. The launch of artificial satellites in the 1960s, which brought new challenges and requirements, has sparked modern analyses of this problem. In complex perturbed environments, the semi-analytical model based on the averaging approach is widely employed as the standard technique in studying secular evolution, providing analytical insights into these perturbing effects, which is quite helpful for mission design. For example, orbital stability analysis (Scheeres et al., 2001; Lara and Juan, 2005; Fu and Wang, 2021), frozen orbit solutions (Delsate et al., 2010; Ulivieri et al., 2013; Condoleo et al., 2016; Circi et al., 2017), and orbits with specific evolution performances (Kamel and Tibbitts, 1973; Liu et al., 2010) have been extensively studied, and relevant results have been applied to a wide range of space missions, as well as lifetime reduction of space debris (Wang and Gurfil, 2016; Wang et al., 2020).

The motion of the hierarchical three-body system, in which the inner two bodies are orbited by an outer body in a wider orbit, is one of the most common problems in astrophysical and astrodynamical contexts. In astrodynamical settings, one of the inner bodies (generally, the spacecraft) can be considered massless, and its motion is dominated by the gravity of the central body and perturbed by the outer body (i.e., the third body or the perturbing body). Considering a highly hierarchical three-body system, the quadrupole level of approximation (also known as the Hill approximation, truncated to the second order of the semimajor axis ratio) is generally retained. In this approximation, based on the framework of the circular restricted-three body problem, both Lidov (1963); Kozai (1962) found that the system is axisymmetric after the double averaging over the inner and outer orbits, leading to the conservation of the ** z**-component of the angular momentum of the inner orbit and, thereby, the integrability of this system. In particular, when the orbit inclination is above the critical inclination (∼

We note that there are other two subcases in which the system exhibits more complex evolutionary behaviors but is still integrable based on a similar mechanism as the Lidov-Kozai ones (i.e., axisymmetry).

1) When the oblateness of the central body is further included but assuming a coplanar geometry (i.e., the perturbing body lies in the central body’s equator) (Kozai, 1969; Scheeres et al., 2001; Delsate et al., 2010; Fu and Wang, 2021). It has been found that the oblateness has changed the phase structures of the orbital dynamics by introducing new equilibrium points, and the critical inclination of the orbital stability changes as a function of the ratio between the perturbing forces of the oblateness and the third-body gravity (namely, the oblateness-related coefficient,

2) The quadrupole-order approximation of the eccentric third-body perturbation (Lidov, 1963; Delsate et al., 2010; Naoz et al., 2013). That is, in the quadrupole-level model, even though the perturbing body is eccentric, the system is integrable. The eccentricity of the perturbing body only induces a small modification in the period of the eccentricity oscillation (De Moraes et al., 2008). However, this is an erroneous conclusion because the quadrupole level approximations have covered up the asymmetry of the system by the exclusion of the argument of perigee of the perturbing body’s orbit (which can be corrected by the octupole order). It has been found that the octupole order can cause a slow, cyclic modulation of the Lidov-Kozai cycles, resulting in orbit flips of the inner orbit (Katz et al., 2011; Lithwick and Naoz, 2011; Naoz, 2016), as well as chaotic motions (Li et al., 2014).

In fact, relaxing either one of the two assumptions (coplanarity or quadrupole-order approximation) will break the integrability of the system and thus lead to qualitatively different evolution behaviors. First, consider an inclined third-body perturbation. Different from the eccentric third-body perturbation problem, there is no special integral of motion in the quadrupole-level model compared to the higher-level models. Actually, only two integrals exist: the semi-major axis and the double-averaged perturbing potential (or the double-averaged Hamiltonian). There are only a few studies contributing to studying the secular evolution including the inclination of the perturbing body, and most of them were based on the quadrupole-order approximation (Yokoyama, 1999; Liu et al., 2012; Circi et al., 2017; Nie et al., 2019; Nie and Gurfil, 2021). Yokoyama (1999) pointed out that the high third-body inclination, with *via* a coordinate transformation. Moreover, since the secular model is based on the quadrupole-order approximation, the critical features of the third-body eccentricity induced by the octupole order are not really included in their models. Then, consider a moderately hierarchical three-body system with an eccentric perturbing body. As stated above, when the perturbing boy’s orbit is eccentric, the quadrupole-order approximation is not enough to describe the secular dynamics and the octupole order is required. In the octupole-level model, the system can be very chaotic and extremely high eccentricities can be excited. In particular, the inner orbit can flip its orientation between prograde and retrograde periodically, accompanied by the eccentricity extremely approaching unity. This secular mechanism is also referred as the eccentric Lidov-Kozai effect (Naoz et al., 2011; Naoz, 2016). Beyond the octupole order, Naoz et al. (2013) point out that higher orders are expected to induce no more dramatic qualitative new features since the integrals of motion remain the same as the octupole order when moving to the higher orders. However, important modifications are possible. Carvalho et al. (2016) extended the perturbing potential to the hexadecapole and dotriocontupole orders and showed the importance of these high-order terms in changing the inversion times of the flip. Will (2017) presented the complete secular equations for both the inner and outer orbits to the hexadecapole order in classical orbit elements. He found that if the system has an equal-mass inner binary, the hexadecapole terms can alone generate orbit flips.

As a consequence, both the inclination and eccentricity of the perturbing body can break the integrability of the classical Lidov-Kozai mechanism. To the best of our knowledge, the two effects have been studied separately, but their joint effects have never been focused. Although several studies considered the eccentricity and inclination of the perturbing body simultaneously, they either focused on the quadrupole level of approximation, actually neglecting the critical effects of the third-body’s eccentricity, or they did not consider the central body’s oblateness, where the third-body’s inclination is actually a problem of coordinate selection.

On the other hand, in the astrophysical context, the discoveries of exoplanetary systems recently have promoted research into this problem (Tokovinin, 2014; Wang et al., 2015; Naoz, 2016; Moe and Kratter, 2021). Compared to the planets in the Solar system with nearly circular and coplanar configurations, the exoplanet systems are found to possess rich geometric architectures, such as eccentric and/or inclined planetary orbits as well as significant stellar oblateness and obliquity (Winn and Fabrycky, 2015). In recent years, dynamical behaviors involving highly eccentric orbits have been invoked to explain a wide range of astrophysical phenomena. For example, at the octupole order (i.e., third order of the semi-major axis ratio), the eccentric outer orbit can drive the inner orbit into extremely high eccentricities, and the sequent tidal dissipation near the pericenter rapidly circularizes that orbit, forming short-period even retrograde planets (Naoz et al., 2011; Naoz et al., 2013; Naoz, 2016). It has also been recently found that the oblateness of a rapidly rotating host star may be responsible for the large mutual inclinations of the exoplanetary systems (Li et al., 2020; Spalding and Millholland, 2020).

Therefore, we see that in a wide array of astrophysical and astrodynamical systems, the two effects can be simultaneously significant and strongly coupled. Concerning this, in this paper, we intend to develop a secular model including the perturbations of the central body’s oblateness and the gravity of an inclined eccentric perturbing body to the hexadecapole order (i.e., the fourth order of the semi-major axis ratio). In particular, the non-spherical gravity of the perturbing body is also incorporated to be applicable for asteroid systems. In order to simplify expressions, we formulate the equations of motion in terms of vectorial elements.

In what follows, we show the applications of our secular model to the binary asteroid system 4951 Iwamoto and a fictitious astrophysical system with assumed parameters to demonstrate its validity and to present some numerical results. Binary asteroid systems are pretty common in the solar asteroid population. About 16 per cent of near-Earth asteroids and main-belt asteroids with diameters below 10 km are binary systems. In recent years, many studies have been devoted to studying the orbital dynamics about these systems, including equilibrium and stability, as well as periodic orbits and invariant manifolds based on the restricted full three-body problem (Bellerose and Scheeres, 2008; Chappaz & Howell, 2015; Dell’Elce et al., 2017; Shi et al., 2019; Shi et al., 2020; Zhang et al., 2020; Li et al., 2021). Recently, Wang and Fu (2020) and Fu and Wang (2021) have extended the semi-analytical orbital theory to an orbiter of a binary asteroid system, assuming a coplanar and circular configuration, by retaining the third-body gravity to the hexadecapole order and including the non-spherical third-body gravity.

This paper is organized as follows. In Section 2, we present the exact Newtonian equations and the double-averaged equations to the hexadecapole order. The process of performing the double averaging, in the form of the vectorial elements, over the inner and outer orbits is elaborated. Next, we extend the hexadecapole-level model to include the non-spherical terms of the third-body gravity in Section 3. Then, we discuss a few applications in Section 4. We compare the different levels of approximation equations with the full equations to validate our semi-analytical model and to show some secular evolutionary behaviors. Finally, we conclude the paper in Section 5.

## 2 Double-averaged model to the hexadecapole order

In this Section, we derive the secular model of the hierarchical restricted three-body system beginning with the full Newtonian equations of motion. By using Legendre polynomials, the disturbing function of the perturbing body is truncated to the hexadecapole order. Then, the average over the inner and outer orbits is performed to filter out the high-frequency variations, arriving at equations that capture the secular orbital evolutions. The averaging results and the final equations are given in terms of the vectorial elements, where the mutual inclination between the central body’s equator and the perturbing body’s orbit is implied in the directions of these vectors.

### 2.1 Full equations of motion

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 *z*-axis along the pole of the central body. For convenient, we represent the perturbing body’s orbital coordinate system as

The perturbation effects of the perturbing body’s gravity on the particle can be stated in the potential form (called the disturbing function) as

where

In hierarchical systems, the disturbing function can be expanded into an infinite series in the form of Legendre polynomials as

where

The disturbing function of the oblateness of the central body can be represented by the gravitational coefficient

where

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

where the disturbing functions

### 2.2 Semi-analytical model

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

where

The validity of the averaging requires that

where

By the averaging theory, the secular Lagrange planetary equations in terms of vector elements, taking the compact and symmetrical form, are given by

where

In the hierarchical three-body systems, as stated above, there are two timescales associated with the orbital motion

where the product of arbitrary two vectors

Substituting Eqs. 11–12 into Eq. 10, disregarding the irrelevant constant terms and scaling to

by virtue of the dyadic identity

Averaging of the octupole-order terms in Eq. 3 requires (Rosengren, 2014)

where the product of three vectors

by dotting

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

Performing the averaging over the disturbing function of the central body’s oblateness, Eq. 4, gives

in which, we note

Then, yields

Finally, by substituting Eqs. 19–23 into the secular Lagrange planetary equations, Eq. 9, the single-averaged equations of motion in terms of

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

The second averaging of the quadrupole-order approximation, Eq. 13, is given as

in which we note

and thus, we have

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

Note that, although we have included the eccentricity of perturbing body’s orbit, there is no components along ** z**-component of the inner orbit’s angular momentum.

The second averaging over the octupole order, Eq. 18, requires

and through which, we have

The secular equations of motion in the octupole order are thus given as

Note that if

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

Finally, applying the averaged disturbing function of the oblateness to the secular equations, we obtain

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

In order to simplify the analysis of the perturbing effects of different perturbing forces, we rewrite the resultant double-averaged disturbing function as

where

To eliminate an extra parameter, we divide the disturbing function by the coefficient

where

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

where

The coefficients

## 3 Secular dynamics with the non-spherical third-body gravity

In this part, we perform an extension of the dynamical model to include the non-spherical gravity of the perturbing body. In the moderately hierarchical three-body systems, such as orbits around one member of a binary asteroid system, the perturbing body’s non-spherical gravity may become important, especially when the eccentricities of the inner and outer orbits bring two of the bodies to close approaches.

### 3.1 The massive two-body model

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))

where

Then, in the synchronous configuration where

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

where

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

### 3.2 Equations with the non-spherical third-body gravity

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

where

in which

For use in perturbation analysis, it is convenient to recast

Substituting

Then, by using Taylor expansion in terms of

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

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

where

## 4 Numerical results

We discuss two typical applications in this Section to demonstrate the obtained secular models and compare their performances at different levels of approximation. First, we apply the secular model to a small-scale binary system, the binary asteroid system, where the system is not highly hierarchical, and the high-order terms and the non-spherical gravity of the perturbing body may be significant. Then, we turn to the large-scale astrophysical systems, where the inclination and eccentricity of the perturbing body may have more critical effects than the hexadecapole terms. By the way, a systematic study of the secular dynamics of such systems is complicated and will be the focus of our future work.

### 4.1 Applications to the binary asteroid systems

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

The second degree and order gravity coefficients in Table 1 are calculated by (Wang and Fu, 2020)

where *s*, *q*, and *p* are the semi-major axes of an ellipsoid, satisfying

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

**FIGURE 3**. Coefficients of *y*-axis), and *y*-axis) with semimajor axis.

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 *via*

**FIGURE 4**. Long-term orbital evolution for **(A–D)** are the temporal evolution of the eccentricity, inclination, longitude of ascending node, and argument of perigee, respectively. The full model is given by Eq. 5 with

**FIGURE 5**. Long-term orbital evolution for **(A–D)** exhibit the temporal evolution of the eccentricity, inclination, longitude of ascending node, and argument of perigee, respectively.

**FIGURE 6**. Long-term orbital evolution for **(A,B)** show the temporal evolution of the eccentricity and inclination, respectively.

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

Remember that in the classical Lidov-Kozai mechanism or the extended “

**FIGURE 7**. Comparisons of the orbital evolution among different combinations of

In Figure 8, we depict the long-term orbital evolution with different inclinations at

**FIGURE 8**. Long-term eccentricity evolution at different inclinations for **(A–D)**, respectively. The integrable case corresponds to the quadrupole-level model with zero third-body inclination (remaining

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

Figure 9 illustrates the case when

**FIGURE 9**. Long-term orbital evolution over 14 scaled time with **(A–C)** show the temporal evolution of the eccentricity, inclination, and argument of perigee, respectively. The evolution of the longitude of node is not shown here. Instead, the variation of the ** z**-component of the orbit angular momentum,

**(D)**.

By using the same initial conditions, we set

**FIGURE 10**. Long-term orbital evolution over 14 scaled time with **(A–D)** show the temporal evolution of the eccentricity, inclination, argument of perigee, and z-component of the orbit angular momentum, respectively.

### 4.2 Applications to other physical systems

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

**FIGURE 11**. Eccentricity and inclination evolution over 100 scaled time for

**FIGURE 12**. Eccentricity and inclination evolution over 100 scaled time for

From Figures 11, 12, we show the important role of the octupole-order terms, especially when the octupole-order coefficient

## 5 Conclusion

Previous studies have revealed the striking features induced by the octupole-order approximation of the third-body perturbation (Katz et al., 2011; Lithwick and Naoz, 2011; Naoz, 2016) and the notable effects caused by the third-body inclination (Yokoyama, 1999; Liu et al., 2012; Nie and Gurfil, 2021), separately, in the hierarchical restricted three-body problem. We find that in some astrophysical and astrodynamical systems, the two factors may be simultaneously significant (Winn and Fabrycky, 2015). Concerning this, in this paper, we have developed a semi-analytical model to provide effective tools for the solution of this problem. Our model includes the perturbations of the central body’s oblateness and the third-body gravity to the hexadecapole order. In particular, the perturbing body is assumed to be in an inclined and eccentric orbit and its non-spherical gravity is also incorporated. Thanks to the double averaging, the secular dynamics can be parameterized by a few coefficients, and through which one can conduct analytical studies on a wide range of physical systems. Formulated in vectorial elements, the equations of motion take a compact form and can be easily reduced to the coplanar or circular perturbing body cases.

Then, numerical integrations have been proposed to demonstrate the validity of our model and to present some numerical results. Comparisons are made among the full, quadruple-level, octupole-level, and hexadecapole-level models. From the numerical results applied to the binary asteroid system 4951 Iwamoto, we have shown the high accuracy of the hexadecapole-level model within the determined valid range of the semimajor axis, while the quadrupole-level and octupole-level models may have obvious deviations. In addition, several secular dynamical features have been highlighted: (a) a joint effect of the perturbing body’s inclination and eccentricity that could drive high eccentricities with chaotic behavior (see Figure 7); (b) the changes in the unstable inclination regions due to the third-body inclination (see Figure 8); and (c) for the

## Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author contributions

TF derived the equations, did the calculations, and wrote the manuscript; YW made the research plan, checked the equations and results, and revised the paper; WH reviewed the paper and gave some suggestions.

## Funding

This work was supported by the National Natural Science Foundation of China (Grant number 11872007).

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2023.1125386/full#supplementary-material

## References

Allan, R. R., and Cook, G. E. (1963). The long-period motion of the plane of a distant circular orbit. *Proc. R. Soc.* 280, 97. doi:10.1098/rspa.1964.0133

Bellerose, J., and Scheeres, D. J. (2008). Restricted full three-body problem: Application to binary system 1999 KW4. *J. Guid. Control. Dyn.* 31, 162–171. doi:10.2514/1.30937

Broucke, R. A. (2003). Long-term third-body effects via double averaging. *J. Guid. Control. Dyn.* 26, 27–32. doi:10.2514/2.5041

Carvalho, J. P. S., Mourão, D. C., de Moraes, R. V., Prado, A. F. B. A., and Winter, O. C. (2016). Exoplanets in binary star systems: On the switch from prograde to retrograde orbits. *Celest. Mech. Dyn. Astron.* 124, 73–96. doi:10.1007/s10569-015-9650-3

Carvalho, J. P. S., Yokoyama, T., and Mourão, D. C. (2022). Single-averaged model for analysis of frozen orbits around planets and moons. *Celest. Mech. Dyn. Astron.* 134, 35. doi:10.1007/s10569-022-10092-6

Chappaz, L., and Howell, K. C. (2015). Exploration of bounded motion near binary systems comprised of small irregular bodies. *Celest. Mech. Dyn. Astron.* 123, 123–149.

Circi, C., Condoleo, E., and Ortore, E. (2017). A vectorial approach to determine frozen orbital conditions. *Celest. Mech. Dyn. Astron.* 128, 361–382. doi:10.1007/s10569-017-9757-9

Condoleo, E., Cinelli, M., Ortore, E., and Circi, C. (2016). Frozen orbits with equatorial perturbing bodies: The case of ganymede, callisto, and titan. *J. Guid. Control. Dyn.* 39, 2264–2272. doi:10.2514/1.g000455

Ćuk, M., and Nesvorny, D. (2010). Orbital evolution of small binary asteroids. *Icarus* 207, 732–743. doi:10.1016/j.icarus.2009.12.005

De Moraes, R. V., Domingos, R. C., and De Almeida Prado, A. F. B. (2008). Third-body perturbation in the case of elliptic orbits for the disturbing body. *Math. Probl. Eng.* 2008, 68.

Dell’Elce, L., Baresi, N., Naidu, S. P., Benner, L. A. M., and Scheeres, D. J. (2017). Numerical investigation of the dynamical environment of 65803 Didymos. *Adv. Sp. Res.* 59, 1304–1320. doi:10.1016/j.asr.2016.12.018

Delsate, N., Robutel, P., Lemaître, A., and Carletti, T. (2010). Frozen orbits at high eccentricity and inclination: Application to mercury orbiter. *Celest. Mech. Dyn. Astron.* 108, 275–300. doi:10.1007/s10569-010-9306-2

Fu, T., and Wang, Y. (2021). Orbital stability around the primary of a binary asteroid system. *J. Guid. Control. Dyn.* 44, 1607–1620. doi:10.2514/1.g005832

Jacobson, S. A. (2011). Dynamics of rotationally fissioned asteroids: Source of observed small asteroid systems. *Icarus* 214, 161–178. doi:10.1016/j.icarus.2011.04.009

Kamel, A., and Tibbitts, R. (1973). Some useful results on initial node locations for near-equatorial circular satellite orbits. *Celest. Mech.* 8, 45–73. doi:10.1007/bf01228389

Katz, B., Dong, S., and Malhotra, R. (2011). Long-term cycling of kozai-lidov cycles: Extreme eccentricities and inclinations excited by a distant eccentric perturber. *Phys. Rev. Lett.* 107, 181101. doi:10.1103/physrevlett.107.181101

Kinoshita, H., and Nakai, H. (1999). *Celest. Mech. Dyn. Astron.* 75, 125–147. doi:10.1023/a:1008321310187

Kinoshita, H., and Nakai, H. (2007). General solution of the Kozai mechanism. *Celest. Mech. Dyn. Astron.* 98, 67–74. doi:10.1007/s10569-007-9069-6

Kinoshita, H., and Nakai, H. (1991). Secular perturbations of fictitious satellites of uranus. *Celest. Mech. Dyn. Astron.* 52, 293–303. doi:10.1007/bf00048489

Kozai, Y. (1969). Long-range variations of orbits with arbitrary inclination and eccentricity. *Vistas Astron* 11, 103–117. doi:10.1016/0083-6656(69)90006-3

Kozai, Y. (1962). Secular perturbations of asteroids with high inclination and eccentricity. *Astron. J.* 67, 591. doi:10.1086/108790

Lara, M., and Juan, J. F. S. (2005). Dynamic behavior of an orbiter around europa. *J. Guid. Control. Dyn.* 28, 291–297. doi:10.2514/1.5686

Lara, M., Rosengren, A. J., and Fantino, E. (2020). Nonsingular recursion formulas for third-body perturbations in mean vectorial elements. *Astron. Astrophys.* 634, A61. doi:10.1051/0004-6361/201937106

Li, G., Dai, F., and Becker, J. (2020). Mutual inclination excitation by stellar oblateness. *Astrophys. J.* 890, L31. doi:10.3847/2041-8213/ab72f4

Li, G., Naoz, S., Holman, M., and Loeb, A. (2014). Chaos in the test particle eccentric kozai–lidov mechanism. *Astrophys. J.* 791, 86. doi:10.1088/0004-637x/791/2/86

Li, X., Qiao, D., and Jia, F. (2021). Investigation of stable regions of spacecraft motion in binary asteroid systems by terminal condition maps. *J. Astronaut. Sci.* 68, 891–915. doi:10.1007/s40295-021-00296-7

Lidov, M. L. (1963). “On the approximated analysis of the orbit evolution of artificial satellites,” in *Dynamics of Satellites/Dynamique des Satellites* (Berlin, Heidelberg: Springer), 168–179.

Lithwick, Y., and Naoz, S. (2011). The eccentric kozai mechanism for a test particle. *Astrophys. J.* 742, 94. doi:10.1088/0004-637X/742/2/94

Liu, X., Baoyin, H., and Ma, X. (2010). Five special types of orbits around mars. *Guid. Control. Dyn.* 33, 1294–1301. doi:10.2514/1.48706

Liu, X., Baoyin, H., and Ma, X. (2012). Long-term perturbations due to a disturbing body in elliptic inclined orbit. *Space Sci.* 339, 295–304. doi:10.1007/s10509-012-1015-8

Lubow, S. H. (2021). An analytic solution to the Kozai–Lidov evolution equations. *Mon. Not. R. Astron. Soc.* 507, 367–373. doi:10.1093/mnras/stab2133

Moe, M., and Kratter, K. M. (2021). Impact of binary stars on planet statistics – I. Planet occurrence rates and trends with stellar mass. *Mon. Not. R. Astron. Soc.* 507, 3593–3611. doi:10.1093/mnras/stab2328

Naoz, S., Farr, W. M., Lithwick, Y., and Rasio, F. A. (2011). Hot Jupiters from secular planet–planet interactions. *Nature* 473, 187–189. doi:10.1038/nature10076

Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., and Teyssandier, J. (2013). Secular dynamics in hierarchical three-body systems. *Mon. Not. R. Astron. Soc.* 431, 2155–2171. doi:10.1093/mnras/stt302

Naoz, S. (2016). The eccentric kozai-lidov effect and its applications. *Annu. Rev. Astron. Astrophys.* 54, 441–489. doi:10.1146/annurev-astro-081915-023315

Nie, T., and Gurfil, P. (2021). Long-term evolution of orbital inclination due to third-body inclination. *Celest. Mech. Dyn. Astron.* 133, 1. doi:10.1007/s10569-020-09997-x

Nie, T., Gurfil, P., and Zhang, S. (2019). Semi-analytical model for third-body perturbations including the inclination and eccentricity of the perturbing body. *Celest. Mech. Dyn. Astron.* 131, 1. doi:10.1007/s10569-019-9905-5

Pravec, P., Scheirich, P., Kusnirak, P., Hornoch, K., Galad, A., Naidu, S. P., et al. (2016). Binary asteroid population. 3. Secondary rotations and elongations. *Icarus* 267, 267. doi:10.1016/j.icarus.2015.12.019

Rosengren, A. J. (2014). Ann and H. J. Smead aerospace engineering Sciences. PhD thesis. (Denver: Univ. Colorado).

Rosengren, A. J., and Scheeres, D. J. (2013). Long-term dynamics of high area-to-mass ratio objects in high-Earth orbit. *Adv. Sp. Res.* 52, 1545–1560. doi:10.1016/j.asr.2013.07.033

Scheeres, D. J., Guman, M. D., and Villac, B. F. (2001). Stability analysis of planetary satellite orbiters: Application to the europa orbiter. *J. Guid. Control. Dyn.* 24, 778–787. doi:10.2514/2.4778

Scheeres, D. J. (2013). Stability of relative equilibria in the full two-body problem. *Ann. N. Y. Acad. Sci.* 1017, 81–94. doi:10.1196/annals.1311.006

Shi, Y., Wang, Y., Misra, A. K., and Xu, S. (2020). Station-keeping for periodic orbits near strongly perturbed binary asteroid systems. *J. Guid. Control. Dyn.* 43, 319–326. doi:10.2514/1.g004638

Shi, Y., Wang, Y., and Xu, S. (2019). Global search for periodic orbits in the irregular gravity field of a binary asteroid system. *Acta Astronaut.* 163, 11–23. doi:10.1016/j.actaastro.2018.10.014

Šidlichovskyý, M. (1983). Quaternions for regularizing celestial mechanics: The right way. *Celest. Mech.* 29, 295.

Spalding, C., and Millholland, S. C. (2020). Stellar oblateness versus distant giants in exciting kepler planet mutual inclinations. *Astron. J.* 160, 105. doi:10.3847/1538-3881/aba629

Tokovinin, A. (2014). From binaries to multiples. II. Hierarchical multiplicity of F and G dwarfs. *Astron. J.* 147, 87. doi:10.1088/0004-6256/147/4/87

Ulivieri, C., Circi, C., Ortore, E., Bunkheila, F., and Todino, F. (2013). Frozen orbital plane solutions for satellites in nearly circular orbit. *J. Guid. Control. Dyn.* 36, 935–945. doi:10.2514/1.59734

Wang, H. S., and Hou, X. Y. (2020). On the secondary’s rotation in a synchronous binary asteroid. *Mon. Not. R. Astron. Soc.* 493, 171–183. doi:10.1093/mnras/staa133

Wang, J., Fischer, D. A., Horch, E. P., and Xie, J. W. (2015). Influence of stellar multiplicity on planet formation. Iii. Adaptive optics imaging of*kepler*stars with gas giant planets. *Astrophys. J.* 806, 248. doi:10.1088/0004-637x/806/2/248

Wang, Y., and Fu, T. (2020). Semi-analytical orbital dynamics around the primary of a binary asteroid system. *Mon. Not. R. Astron. Soc.* 495, 3307–3322. doi:10.1093/mnras/staa1229

Wang, Y., and Gurfil, P. (2016). Dynamical modeling and lifetime analysis of geostationary transfer orbits. *Acta Astronaut.* 128, 262–276. doi:10.1016/j.actaastro.2016.06.050

Wang, Y., Luo, X., and Wu, X. (2020). Long-term evolution and lifetime analysis of geostationary transfer orbits with solar radiation pressure. *Acta Astronaut.* 175, 405–420. doi:10.1016/j.actaastro.2020.06.007

Will, C. M. (2017). Orbital flips in hierarchical triple systems: Relativistic effects and third-body effects to hexadecapole order. *Phys. Rev. D.* 96, 023017. doi:10.1103/PhysRevD.96.023017

Winn, J. N., and Fabrycky, D. C. (2015). The occurrence and architecture of exoplanetary systems. *Annu. Rev. Astron. Astrophys.* 53, 409–447. doi:10.1146/annurev-astro-082214-122246

Yokoyama, T. (1999). Dynamics of some fictitious satellites of Venus and Mars. *Space Sci.* 47, 619–627. doi:10.1016/s0032-0633(98)00110-x

Keywords: semi-analytical model, inclined eccentric third-body perturbation, eccentric Lidov-kozai mechanism, non-spherical third-body gravity, hexadecapole order

Citation: Fu T, Wang Y and Hu W (2023) Semi-analytical orbital model around an oblate body with an inclined eccentric perturber. *Front. Astron. Space Sci.* 10:1125386. doi: 10.3389/fspas.2023.1125386

Received: 16 December 2022; Accepted: 07 February 2023;

Published: 23 February 2023.

Edited by:

Ludmila Carone, Austrian Academy of Sciences, AustriaReviewed by:

Yu Jiang, China Xi’an Satellite Control Center, ChinaJean Paulo Carvalho, Federal University of the Recôncavo of Bahia, Brazil

Copyright © 2023 Fu, Wang and Hu . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Yue Wang, ywang@buaa.edu.cn