An Orbit-flip Mechanism by Eccentric Lidov–Kozai Effect with Stellar Oblateness

Orbit flips have been previously found under the eccentric Lidov–Kozai effect (ELK) in hierarchical three-body systems. Recently, we have found that, in certain conditions, the orbit can flip its orientation in a much different manner, where the stellar oblateness plays an important role. In this paper, orbit-flip behaviors with the ELK effect are investigated as the stellar oblateness varies within a wide range. This is of significance because recent works have shown that the oblateness of young stars has a widespread distribution and may have critical effects on sculpting the final orbital states of close-in planets. Our dynamical model includes the secular potential of the perturber to octupole order and the secular effects of the stellar oblateness. An alignment between the orbit plane of the outer perturber and the stellar bulge is assumed. Our findings mainly consist of two aspects. (i) A new type of orbit-flipping mechanism induced by a combination of the ELK and stellar oblateness effects, referred to as the ELK–J 2 effect, is discovered and confirmed. (ii) We demonstrate that, in the considered aligned configuration, the stellar oblateness suppresses orbit flips due to the ELK effect and produces new flips through the ELK–J 2 effect. Moreover, if the stellar oblateness perturbations are of the same order as the octupole perturbations of the outer perturber, the ELK-induced orbit flips are almost entirely suppressed, while the ELK–J 2 effect reaches its peak for the considered strength of the octupole perturbations. However, from a global view, stellar oblateness always reduces flipping orbits.


Introduction
Three-body systems are common in astrophysics, from planet-satellite systems to triple-star and black hole systems, covering a wide range of scales and configurations.Most of them are hierarchical for reasons of stability, i.e., the inner binaries are orbited by an exterior companion on a much wider orbit (Naoz 2016).Although the general three-body problem is mathematically nonintegrable, the motion of the hierarchical triples can be described by two slowly evolving Keplerian orbits under their secular interactions, and then various analytical approximations can be applied.Lidov (1962) and Kozai (1962) studied the hierarchical threebody problem by retaining the expansion of the interaction potential to the quadrupole order.They considered applications when one of the inner bodies is massless, i.e., when the outer orbit contains most of the angular momentum of the system.After double-averaging over the inner and outer orbits, the secular potential is axisymmetric, and the longitude of the ascending node becomes a cyclic variable in the Hamiltonian model, indicating an additional constant of motion, the projected angular momentum (j e i 1 cos z 2

= -
).As a consequence, the secular dynamics are integrable and can be solved analytically (Kinoshita & Nakai 1999, 2007).Notably, it is found that when the mutual inclination between the inner and outer orbits is beyond the critical inclination (39.2°), the eccentricity and inclination undergo largeamplitude, periodic oscillations on a secular timescale, known as Lidov-Kozai cycles (LK cycles).The timescale of the LK cycles, t LK , at the quadrupole order is approximated by (Antognini 2015)  where m 0 and m per are the masses of the central object and the outer perturber, P in and P out are the periods of the inner and outer orbits, and e per is the eccentricity of the perturber.In recent years, the LK mechanism has been widely invoked to study the formation and evolution of a wide variety of astrophysical systems (Eggleton & Kiseleva-Eggleton 2001;Wu & Murray 2003;Fabrycky & Tremaine 2007;Perets & Fabrycky 2009;Antonini et al. 2014;Dong et al. 2014;Naoz & Fabrycky 2014;Storch et al. 2014;Shevchenko 2017;Bhaskar et al. 2021).
It was recently recognized that, if the octupole order is included for an eccentric outer orbit, a system with an eccentric Lidov-Kozai effect (ELK) exhibits new and striking features that are impossible in the standard LK mechanism, including arbitrarily high eccentricities, chaos, and orbit flips (Katz et al. 2011;Lithwick & Naoz 2011;Naoz et al. 2011Naoz et al. , 2013;;Li et al. 2014a).The octupole order of the outer perturber's potential causes the LK cycles to undergo slow evolution on a very-longterm timescale, t ELK , where t ELK ?t LK (Katz et al. 2011;Lithwick & Naoz 2011).This happens because the octupole order breaks the axisymmetry of the secular dynamics, and then the projected angular momentum is no longer constant.According to Antognini (2015), t t , where ε oct measures the importance of the octupole order relative to the quadrupole order.
One of the most noteworthy aspects of the ELK effect is that the inner orbit can flip its orientation between prograde and retrograde, with the eccentricity arbitrarily approaching unity.In the test-particle limit, Katz et al. (2011) averaged the doubleaveraged secular equations of motion over the LK cycles for the first time and identified a new constant of motion.By applying this new constant, they established an analytical criterion for orbit flips.For the same problem, Lithwick & Naoz (2011) used a numerical technique to investigate the dynamical structure and flipping regions of the ELK effect.Li et al. (2014b) found that, in addition to the orbit flips with inclinations within the range 40 , 140    [ ] , characterized by initial low eccentricities and high inclinations, there is another type of flip that can occur in a near-coplanar configuration (i 0 ≈ 0°or 180°), characterized by initial high eccentricities and low inclinations.Lei (2022a) confirmed that the above two types of flips are dynamically similar but located in different regions of the phase space.In addition, there exists a third type, in which orbit flips happen with initial intermediate eccentricities.Sidorenko (2018) and Lei & Gong (2022b) demonstrated that the ELK effect is essentially a secular resonance and the flipping orbits are resonant trajectories with the resonant center located at i = 90°.
Evidence shows that weak short-range effects, which are significant during the close pericenter passages of the LK cycles, may play an important role in forming the final orbits of exoplanets in the process of high-eccentricity migration.From a geometric view, the growth of eccentricity depends on the position of the apsidal line (Fu & Wang 2021).Therefore, the short-range effects that induce the pericenter precession, including the general relativity, rotational bulge, and tidal distortion, can significantly suppress the LK cycles (Wu & Murray 2003;Fabrycky & Tremaine 2007;Dong et al. 2014).Furthermore, these short-range effects could impose an impact on orbit flips and limit the maximum achievable eccentricity induced by the ELK effects (Liu et al. 2015).
On the other hand, recent works have demonstrated that the oblateness or rotational bulge of young stars could have a comparable or even dominant influence on the close-in planets relative to the exterior perturber (Spalding & Batygin 2016;Li et al. 2020;Spalding & Millholland 2020).The oblateness, J 2 , of a star deformed by its own rotation can be estimated by (Sterne 1939;Spalding & Batygin 2016) where a W • is the stellar rotational frequency, Ω b is the break-up rotational frequency, and k 2 is the classical apsidal motion constant.As the star spins down, J 2 decays linearly with time over a very long timescale.This implies that the stellar oblateness may be an important factor in forming close-in planets and hot Jupiters at the early stage of the system (Spalding & Batygin 2016;Li et al. 2020).Although the stellar oblateness has been included as one of the weak short-range forces in some recent works, a systematic study of the LK/ ELK effect combined with J 2 varying within a wide range is currently lacking.
In this paper, the secular evolution of a test particle orbiting around an oblate central star with an eccentric exterior companion is studied.That is, two perturbing effects are included: the stellar oblateness and the tidal gravitational force from the outer companion.The secular interaction potential between the inner and outer orbits is expanded up to octupole order, and the stellar oblateness is chosen within a wide range.
The efforts are concentrated on a new orbit-flip mechanism with this dynamical context, i.e., the ELK-J 2 effect.By combining analytical and numerical investigations, the underlying mechanism, global phase space structures, and dynamical features during the secular evolution are systematically studied.
The rest of the paper is organized as follows.Section 2 presents the secular equations of motion in a vectorial form.In Section 3, the new orbit-flip mechanism with the ELK-J 2 effect and its dynamical essence are introduced as an extension of the LK-J 2 effect.In Section 4, by employing Poincaré sections, the global structure of the phase space is uncovered, and the flipping regions are determined therein.Finally, the paper is concluded in Section 5.

Equations of Motion
Consider a massless test particle (an inner planet) orbiting around an oblate central star of mass m 0 in the presence of a massive exterior perturber of mass m Per .In the test-particle limit, the perturber is moving on a fixed eccentric Keplerian orbit coplanar with the stellar equatorial plane, and the particle is moving on a slowly evolving orbit perturbed by the stellar oblateness and the gravity of the perturber.Note that the secular effects of the inner orbit on the stellar rotation and the outer orbit are ignored due to the massless assumption.The explicit discussions about this assumption are elaborated in Appendix B. Specifically, for a solar-mass star with rotationinduced oblateness J 2 ≈ 10 −4 and with a planet orbiting at 1 au, the planet's mass m should have m  10m Earth to disregard its effect on the stellar spin.
It is convenient to describe the particle's orbit by two vectors: the dimensionless angular momentum j j e where j ˆis the unit vector along j, and the eccentricity e e e = ˆ, where the unit vector e ˆis pointing in the direction of the pericenter.They satisfy j • e = 0 and j 2 + e 2 = 1.In the same way, j Per and e Per are the geometric representations of the perturber's orbit.The geometry of the hierarchical three-body system is shown in Figure 1.
The inertial coordinate system is defined by selecting the perturber's orbital plane as the invariant plane, where the origin is located at the central star, the x-axis and z-axis are aligned with e Per and j Per , respectively, and the y-axis completes the orthogonal right-handed reference frame.In the coplanar configuration, the stellar spin axis n 0 ˆlies along the z-axis.To capture the secular evolution of the particle's orbit, averaging over the periods of the inner and outer orbits is required.Following a standard averaging procedure, the single averaging of the stellar oblateness potential over the inner orbit is given by (Wang & Fu 2020) As for the perturber's potential, the double-averaged quadrupole and octupole-order terms are given by (Liu et where J 2 and R 0 are the oblateness and radius of the central star, and κ Ob and ε Oct measure the stellar oblateness perturbations and the outer perturber's octupole perturbations relative to the outer perturber's quadrupole perturbations, respectively.By using these averaged potentials, the secular Lagrangian equations of motion, in terms of the vectorial elements j and e, take an elegant and symmetrical form as (Allan & Cook 1964)  The explicit expressions of the secular equations of motion are given in Appendix A.

A New Type of Orbit Flip with the ELK-J 2 Effect
The striking phenomena of orbit flips was first discovered in an octupole-level approximation with an eccentric outer orbit.The related mechanism is referred as the ELK effect (Lithwick & Naoz 2011;Naoz 2016).Recent research has demonstrated that the ELK effect is essentially a secular resonance with the critical argument j sign z v w = W + ( ) , i.e., the longitude of pericenter for prograde or retrograde orbits, which librates around 0 or π on a very long timescale (Sidorenko 2018;Lei 2022a;Lei & Gong 2022b).In this section, we are going to reveal a new orbit-flip mechanism by the combined ELK and stellar oblateness (J 2 ) effects.The exploration of the dynamical essence will be presented as an extension of the LK-J 2 effect studied before.

An Example of an Orbit Flip by the ELK-J 2 Effect
The left column of Figure 2 illustrates an example of an orbit flip by the ELK-J 2 effect.For comparison, an example of an orbit flip by the ELK effect in the absence of stellar oblateness is shown in the right column.The system considered is characterized by κ Ob ≈ 0.1775 and ε Oct = 0.0393 for the ELK-J 2 case, and κ Ob = 0 and ε Oct = 0.0393 for the ELK case.
In the left column of Figure 2, i.e., under the ELK-J 2 effect, the orbit periodically flips its orientation, the inclination changes in an interval of about 70 , 110   ( ) centered at i = 90°, the eccentricity oscillates between 0.67 and 0.88, and both the arguments ω and Ω librate around 0. Since Ω changes on a much longer timescale and with a larger amplitude than ω, the libration of the arguments ω ± Ω has a similar timescale as Ω.
In contrast, the system under the ELK effect in the right column behaves in a different manner.The eccentricity is excited to extremely high values with e 1 max on the order of 10 −5 and the inclination can oscillate quite far away from 90°.Both the arguments ω and Ω circulate within 0, 2p ( ), although Ω changes its circulating direction at some specific points.The comparable frequencies of ω and Ω lead to a quite long timescale of the libration of the argument j sign z v w = W + ( ) .In conclusion, the two types of orbit flip presented in Figure 2 are characterized by different types of ω motion: ω is librating around ω = 0 in the ELK-J 2 case, while it is circulating in the ELK case.Therefore, an orbit flip by the ELK-J 2 effect should have a qualitatively different mechanism from that by the ELK effect.Note that the mode of ω motion is mainly determined by the quadrupole-level approximation.Under the standard LK mechanism at the quadrupole level, ω is librating around ω = π/2 (or 3π/2) or circulating, but there are no such librating cycles around ω = 0 (or π).Therefore, the new type of ω motion in the ELK-J 2 effect is produced by further including the stellar oblateness in the LK mechanism (San-Juan et al. 2006;Delsate et al. 2010;Fu & Wang 2021).At the quadrupole-level approximation, this dynamical effect is referred to as the LK-J 2 mechanism, since the eccentricity of the outer orbit would not show up.

The ELK-J 2 Mechanism as a Secular Resonance
At the quadrupole-level approximation, j z is conserved in the considered aligned configuration.This constant gives rise to the integrability of the system.As a result, the cyclic motion of the system can be fully described by two constants: j z and the energy constant ψ ≈ <ψ Ob > + <ψ Quad >.
Figure 3 depicts nontrivial phase portraits of a system with κ Ob = 0.1775 by contour plots of the energy constant for the given j z .The phase space structure shown in panel (A) is inherited from the standard LK mechanism.The two modes of ω motion, circulation and libration, are known as LK cycles.When j z approaches 0, i.e., i comes close to 90°, the motion will transition to the phase space (B) and (C), where new structures appear.Librating cycles around ω = 0 and π and circulating cycles from ω = 2π to ω = 0 emerge, which are referred to as LK-J 2 cycles.Both the LK and LK-J 2 cycles are so-called quadrupole-level cycles.
When moving to an octupole-level approximation, if the perturber is eccentric there is only one constant of motion left, i.e., the secular energy ψ.Therefore, the secular evolution of the particle's orbit is more complicated and even chaotic.The octupole order can force the quadrupolelevel cycles to undergo very-long-term modulations at a timescale of t ELK .On a shorter timescale of t LK , the system evolution can still be described by the quadrupole-level cycles.Therefore, the system behavior in the phase space can be determined if we could capture the very-long-term evolution of j z .This process is similar to the investigations of the ELK effect by Lithwick & Naoz (2011) and Katz et al. (2011).Bearing this idea in mind, we present a qualitative analysis of the dynamical essence of the ELK-J 2 effect below.
For the example shown in the left column of Figure 2, the system is initialized on a librating LK-J 2 cycle around ω = 0 in the phase space (B), indicated by a magenta star in Figure 3.Over the integration time, ω librates around 0, and Ω also librates around 0 but on a much longer timescale.The specific evolutionary behaviors of ω and Ω will be shown to be relevant to orbit flips induced by the ELK-J 2 effect.
Let us take a closer look at this issue.Since the system is in quadrupole-level cycles with j z = 1, by taking the lowest-order terms of j z , the change of j z derived from Equation (A9) is The equation for Ω, at the quadrupole-dominated level, is Note that  W changes sign when j z crosses zero at the orbit flip.From Equation (12), if j does not come close to 0, i.e., e is adequately smaller than 1, Ω changes slowly due to the small value of j z .On the other hand, for the librating cycles close to ω = 0 and ω = π, ω is librating with small amplitudes.If the amplitude of ω is much smaller than that of Ω, both the arguments Ω + ω and Ω − ω in Equation (11) will vary slowly.Therefore, according to Equation (11), j z  also changes slowly.
Then, it may provide enough time for j z to cross zero before j z  changes its sign, i.e., before the orbit flips its orientation, as shown in Figure 2.Under this mechanism, orbit flips are based on the librating LK-J 2 cycles around ω = 0 and ω = π, which provide small amplitudes of ω and a limited maximum eccentricity when i is close to 90°.A numerical confirmation of the above analyses is carried out in the e , 0 0 w ( ) space in Figure 4, where each dot in the space denotes a unique trajectory.Trajectories that flip their orientations during the evolution are marked by red.Considering the phase portraits of the quadrupole-level approximation when j z is small in Figure 3, it is clear that the flipping orbits are limited in a region near the libration center ω = 0, which is consistent with our above analyses.
As for the orbital flips due to the ELK effect, in the case of the standard LK cycles (κ Ob = 0) with j z = 1 the eccentricity is excited to extremely high values, i.e., j → 0, during each cycle.According to Equation (12), Ω changes fast in this case.For the circulating LK cycles, if the frequencies of ω and Ω are close to each other, their changes over time can cancel each other out under certain conditions, resulting in a slow variation of j sign z w W + ( ) , as shown by the right column of Figure 2.This secular resonance has been confirmed to be the underlying mechanism of the ELK effect (Sidorenko 2018; Lei & Gong 2022b).
From the above, we can conclude that the ELK-J 2 effect is a secular resonance with critical arguments of ω ± Ω, based on the librating LK-J 2 cycles around ω = 0 or π, whereas the ELK effect is a secular resonance with the critical argument of j sign z v w = W + ( ) , based on the circulating LK cycles from ω = 0 to ω = 2π.Since these two dynamical effects occur in different regions of the phase space, the distinct quadrupolelevel cycles cause qualitative differences between the ELK and ELK-J 2 effects, and the orbital flips induced by them.
One of the notable differences is that the orbit flips induced by these two mechanisms have different inclination and eccentricity coverages.The secular evolutions of j z , e, and i for the two examples in Figure 2 are shown in Figure 5.As we can see, although the two trajectories evolve within similar intervals of j z , they actually evolve in different regions of the e i , ( ) space.The LK-J 2 cycles (left column) are restricted to lower eccentricities and, thus, a narrow inclination interval  w ( ) space for the system presented in Figure 2 (κ Ob ≈ 0.1775, ε Oct ≈ 0.0393) by numerical simulations.For each grid, the initial Ω is Ω 0 = 0 and the initial i is given by j z,0 = 0.105 and 0.05 for the left and right panels, respectively.The initial conditions that lead to orbit flips are marked by red.centered at 90°, whereas the LK cycles (right column) are excited to extremely high eccentricities accompanied by an extensive inclination range when e → 1.So far, a preliminary qualitative analysis of the dynamical essence of the ELK-J 2 effect and the orbital flips induced has been conducted.Next, a numerical method will be adopted to investigate the global structure of the phase space of the ELK-J 2 effect and to determine the orbit-flipping regions in the e i , 0 0 ( ) space.

Global Structure of the Phase Space
In this section, the Poincaré surface of section is exploited to explore the global structure of the phase space of the ELK-J 2 effect.The Poincaré section is a powerful tool for studying dynamical systems with a high-dimensional phase space, such as the ELK effect in hierarchical triple systems (Lithwick & Naoz 2011;Li et al. 2014a;Naoz et al. 2017;Lei 2022a).The section surfaces at different energy levels can capture all the phase space structures.Through these sections, the orbitflipping boundaries in the parameter space can be determined by the boundaries of the libration islands centered at 90°.
As highlighted in Section 3, the ELK-J 2 effect is based on the librating cycles around ω = 0 or π.Accordingly, the Poincaré section is defined by 0, 0 and 0, 0. 13 Thanks to the symmetry, the librating cycles at ω = π have symmetrical structures with those at ω = 0. Therefore, only the results about librating cycles at ω = 0 will be presented. ),are chosen to present the phase space structure.The range of the energy constant in the section ω = 0 is given in Figure 6, where only the energy curves crossing i = 90°are considered.The Poincaré sections are reported in Figure 7.
In Figure 7, the Poincaré sections of the upper panels are defined by 0, 0  w w = > , and in the lower panels they are defined by 0, 0  w w = < .The libration islands centered at i = 90°dominate the phase space structure of the ELK-J 2 effect.In each panel, the green dots indicate the boundaries of these libation islands within which orbit flips will occur.
The phase space structure changes with the values of the energy constant.In panel (a), with ψ * = 0.42, there is a single libration island centered at i = π/2, Ω = π and a single saddle point located at i = π/2, Ω = 0.The curve through the saddle point separates the flipping and nonflipping regions.
When the energy constant decreases to ψ * = 0.34, as shown in panel (b), the libration center at Ω = π bifurcates into two libration centers and a saddle point, all located at i = 90°.As the energy constant continues to diminish, the two libration centers move toward Ω = 0, and the trajectory through the saddle point becomes a new dynamical separatrix.Finally, it changes into the structure shown in panel (c).
The sections in panels (d) and (e) are mapped with the same energy levels as in panels (c) and (a), respectively.Indeed, each pair with the same energy constant stands for the same group of trajectories.As shown in the lower panels of Figure 3, the librating LK-J 2 cycles cross the ω = 0 sections twice within one period, one with 0  w < and the other with 0  w > .Since the librating LK-J 2 cycles in the e , w ( ) plane are counterclockwise, 0  w > sections give the minimum eccentricity and 0  w < capture the maximum eccentricity during each cycle.
Panel (f) represents the phase space structure at extremely high eccentricities, i.e., the regions with high energy levels in Figure 6.In the current system, for orbits with extremely high eccentricities, the strong effect of the stellar oblateness during the close pericenter forces them to undergo trivial evolutions.Therefore, we would not pay much attention to the orbits in this region.
By analyzing the libration islands at different levels of energy constant, the flipping regions in the entire phase space can be determined.Lei (2022a) recently adopted a similar strategy to study the flipping regions under the ELK effect.In Figure 7, i cos D at the libration center is utilized as the width of the libration island to represent the boundaries of flipping regions in the e i , 0 0 ( ) space.The corresponding e 0 at these boundaries can be obtained by numerical integrations or by exploiting the energy constant.
The results are reported in Figure 8, where several distinct flipping regions can be observed; the dominant two are designated as I and II, respectively.Region I, the low-e 0 area, denotes the libration islands centered at Ω = π, corresponding to the phase space structure shown in panels (a) and (b) of Figure 7.In this region, orbit flips are confined to a narrow interval symmetric about i = 90°.Region II is given by the libration islands centered at Ω = 0, i.e., panels (c) and (d) of Figure 7.This region dominates the ELK-J 2 effect, allowing a relatively wide range of flipping e 0 and i 0 .The other bounded regions of Figure 8, occupying the high-e 0 areas, determined by the phase space structure similar to panels (e) and (f), are trivial due to their small size in e 0 or i 0 .
Then, an exploration of flipping regions for orbits initialized in the entire e i , 0 0 ( ) space with given initial ω 0 and Ω 0 is performed by numerically integrating Equations (A5)-(A10).All the trajectories that undergo orbit flips over the integration time, 500t LK , are recorded in Figure 9.It indicates that the numerical results agree well with the analytical results in Figure 8.However, in the range of e 0.7, 0.95 0 = ( ) , the numerical integrations define a slightly larger inclination ), in the e i , cos ( ) plane with ω = 0 for a system with κ Ob = 0.1.For clarity, we adopt a different abscissa scale for the e > 0.9 areas.interval than the analytical approach.This is attributed to the chaotic orbit-flipping behaviors captured by the numerical method but not included in the analytical results.Further, chaotic motions enable orbit flips at Ω 0 = π (red dots) in the regions dominated by the libration islands centered at Ω 0 = 0 (black dots).

Phase Space with Different Values of κ Ob
Numerical explorations are here extended to involve various κ Ob .First, with ε Oct = 0.03, κ Ob is assumed to be 0, 10 −4 , 10 −3 , and 0.03 to examine the changes of the phase space characterized by different strengths of the stellar oblateness.Among these, κ Ob = 0 represents the ELK effect in the absence of stellar oblateness.The results are reported in Figure 10.
The ELK effect and ELK-J 2 effect can be distinguished by the modes of ω motion, i.e., the quadrupole-level cycles.The former is characterized by the circulation of ω from 0 to 2π, i.e., circulating LK cycles, whereas the latter is characterized by the libration of ω around 0 or π, i.e., librating LK-J 2 cycles.As shown in Figure 10, the two effects can simultaneously appear in the phase space of a system with appropriate κ Ob and ε Oct .
In the case of the ELK effect with κ Ob = 0 and ε Oct = 0.03, in the upper-left panel of Figure 10  space with low, intermediate, and high e 0 , respectively.As the phase space structure of the ELK effect is qualitatively different from the ELK-J 2 effect shown in Figure 7, there is no connection between the flipping regions shown in this panel and that shown in Figure 8.One of the interesting features of the ELK effect is that in the high-e 0 region, when e 0 is larger than a critical value orbit flips can occur from a nearly coplanar configuration (i.e., i 0 ∼ 0°and 180°).This phenomenon was first discovered by Li et al. (2014b).
In the case of κ Ob = 10 −4 , in the upper-right panel of Figure 10, although the stellar oblateness is quite weak, the structures of the phase space change a lot.Notably, the ELK-J 2 effect arises near i = 90°, while the ELK effect near i = 90°is suppressed.In the intermediate-e 0 parts, there is an inclination gap between the ELK-and ELK-J 2 -dominated regions, where orbit flips are completely suppressed over the total integration time.Indeed, the gap also exists in the high-e 0 region, although it is tiny in this case.Furthermore, the near-coplanar orbit flips are now impossible, and the flipping inclination interval shrinks when e 0 → 1.
In the case of κ Ob = 10 −3 , in the lower-left panel of Figure 10, the gaps become larger.The ELK effect is further suppressed and the ELK-J 2 effect grows more significant.
Finally, when κ Ob = 0.03, in the lower-right panel of Figure 10, the ELK effect almost vanishes in the entire e i , 0 0 ( ) space and the ELK-J 2 effect dominates the structure of the phase space.
Figure 10 shows how the structure of the entire e i , 0 0 ( ) space changes as the strength of the stellar oblateness increases.Next, we examine these changes in more detail, examining the suppression effect of stellar oblateness on the ELK effect and the development of the ELK-J 2 effect.
In Figure 11, each curve marks the minimum flipping inclination, i 0 , for the given ε Oct and different values of κ Ob .Below these curves (lower i 0 ), there are no flipping orbits, but above them, as i 0 increases, a growing number of orbits flip.
The curve in Figure 11 rises as the strength of the stellar oblateness (κ Ob ) grows, until it reaches a turning point.The turning point is attributed to the transition of the phase space to the ELK-J 2 effect.At this point, the ELK effect is almost completely suppressed.The turning points of the curves in Figure 11 with ε Oct = 0.01, 0.03, 0.06, and 0.1 are given by κ Ob ∼ 0.003, 0.032, 0.068, and 0.136, respectively.Therefore, if ε Oct lies within (0.01, 0.1), we can arrive at a rough estimate on the suppression effect of stellar oblateness on the ELK effect: when the stellar oblateness is of the same order as the octupole terms, i.e., κ Ob ∼ ε Oct , orbit flips by the ELK effect are almost completely suppressed.
Analytically, recall that orbit flips due to the ELK effect rely on the secular resonance with the critical argument j sign z w W + ( ) .However, the pericenter and nodal precession induced by the stellar oblateness would affect this resonance; that is, the strength of stellar oblateness measured by κ Ob matters.On the other hand, the timescale of the evolution of j z depends on the strength of the octupole terms, or, equivalently, the timescale of the ELK effect, characterized by ε Oct , as indicated by Equation (11).Therefore, it is theoretically possible to derive a relation between κ Ob and ε Oct to describe such suppression effects.
This conclusion (κ Ob ∼ ε Oct ) is achieved by numerical integrations with the initial e 0 close to zero, but it can be used to approximately represent all initial e 0 cases, as is confirmed by Figure 10. ( ) space for systems with κ Ob = 0, 10 −4 , 10 −3 , and 0.03, respectively.The initial conditions are indicated in these plots.The total integration time is 500t LK for each trajectory.
Figure 12 presents the generation of orbit flips by the ELK-J 2 effect.We follow a similar process as in Figure 11 to gain these curves but with initial e 0 = 0.82 and Ω 0 = 0 (ω 0 = 0 remains unchanged).For a given e 0 and Ω 0 , orbit flips occur in the region dominated by the ELK-J 2 effect.However, note that the curves in this plot do not represent the minimum flipping i 0 for all the values of Ω 0 (see the caption for explicit description).
It is interesting to note that the curves of different ε Oct in Figure 12 follow a similar regular decreasing curve at the beginning, and then as κ Ob increases a separation point suddenly arises for each curve.The separation points appear sequentially in the order of the values of ε Oct .From this plot, for ε Oct = 0.01, 0.03, 0.06, and 0.1, the separation points are given by κ Ob ∼ 0.004, 0.039, 0.069, and 0.10, respectively.For all the values of κ Ob , each curve reaches its minimum at the separation point for the given ε Oct .Therefore, roughly, we conclude that when the stellar oblateness perturbations are of the same order as the octupole perturbations, i.e., κ Ob ∼ ε Oct , the ELK-J 2 effect becomes most significant for the given ε Oct .
After that, the four curves in Figure 12 finally end at i 0 = 90°( given they are symmetric about i 0 = 90°) with a similar value of κ Ob (∼0.5).This means that, for systems with ε Oct  0.1, almost all flipping orbits would be suppressed when κ Ob  0.5.
In sum, two observations can be made.(i) With the increase of stellar oblateness, the ELK effect is gradually suppressed starting from i 0 = 90°, while the ELK-J 2 effect arises in turn.Moreover, when the stellar oblateness is of the similar order as the octupole terms, i.e., κ Ob ∼ ε Oct , the ELK effect is nearly entirely suppressed; further, in these cases, the ELK-J 2 effect becomes most significant, as shown by Figures 11 and 12. (ii) Near-coplanar orbit flips can be prohibited by a weak stellar oblateness.

Misaligned Configuration
In the above, we have considered systems with alignment between the stellar spin and the angular momentum of the outer orbit.Here, we present a brief examination of a misaligned configuration.In an aligned system, the secular potential is symmetric about the invariant plane, leading to symmetric evolution between the prograde and retrograde orbits.Misalignment breaks such symmetry.Therefore, the secular dynamics is expected to be more complicated and chaotic.
Figure 13 shows a reproduction of the last panel of Figure 10 for misaligned cases, where the stellar obliquities are ϑ = 10°( upper panel) and 40°(lower panel), respectively.From this picture, we can observe that (i) when the stellar obliquity is small (ϑ < 10°), the system exhibits similar flipping behaviors as the aligned one; (ii) when the stellar obliquity is significant (ϑ ∼ 40°), the structure changes a lot; further, the ELK and the ELK-J 2 effects in the phase space are almost indistinguishable; lastly, (iii) it is interesting that the stellar oblateness does not always suppress orbit flips.As in the case of ε Oct = 0.03, the minimum flipping i 0 is ∼50°when e 0 ∼ 0 for κ Ob = 0.03 and ϑ = 40°(the lower panel of Figure 13), while the value is ∼60°f or κ Ob = ϑ = 0 (the upper-left panel of Figure 10).That is, in this case the stellar oblateness widens the flipping inclinations.9) and (10) initialized with e 0 = 0.01, ω 0 = 0, and Ω 0 = π, scanning through the values of i 0 from 0 to π/2 to determine the minimum i 0 that allows flips.We did not scan the initial Ω 0 , because according to the phase space structure at e 0 ∼ 0, the lower limit of flipping i 0 is determined at Ω 0 = π.At this initial e 0 , orbit flips occur in the region dominated by the ELK effect.These curves are symmetric about i 0 = 90°, and only the parts below 90°are depicted.
Figure 12.Minimum flipping inclination due to the ELK-J 2 effect.These curves are drawn following a similar approach as in Figure 11, but setting initial e 0 = 0.82 and Ω 0 = 0 (ω 0 = 0).This picture captures orbit flips mostly due to the ELK-J 2 effect, as indicated by Figure 10.However, note that the curves in this figure do not represent the minimum flipping i 0 for all initial Ω 0 , because at low values of κ Ob flips due to the ELK effect with libration islands centered at Ω 0 = π generally determine a lower minimum flipping i 0 .However, at high values of κ Ob , when the ELK effect has been almost entirely suppressed, the initial Ω 0 = 0 gives the exact minimum values.
We highlight that the above conclusions are not comprehensive.A systematical study of the misaligned configuration is complicated and deserves future investigation.

Conclusions
In this paper, secular orbital evolution is investigated under the combined effects of stellar oblateness and an eccentric distant perturber in a hierarchical three-body system.The perturbations of the outer perturber is formulated up to the octupole order.At this octupole-level approximation, it is found that the stellar oblateness has a suppression effect on orbit flips due to the ELK effect.In addition, it meanwhile gives rise to a new kind of orbit flip, which we refer to as the ELK-J 2 effect.The main results are summarized as follows: 1.A new type of orbit flip induced by the ELK-J 2 effect arises when the stellar oblateness is significant.The ELK-J 2 effect is found to be a secular resonance with critical arguments of ω ± Ω.In contrast, the ELK effect has been proven to be a secular resonance with the critical argument of j sign z w W + ( ) .Essentially, they are based on different quadrupole-level cycles: the former is attributed to librating LK-J 2 cycles around ω = 0 or π, and the latter is the result of circulating LK cycles with ω from 0 to 2π. 2. A global map of the phase space structure reveals that the libration islands centered at i = π/2, which allow orbit flips, dominate the phase space structure of the ELK-J 2 effect.By stepping through energy levels, the flipping boundaries in the e i , 0 0 ( ) space can be determined.3. Increases of the strength of the stellar oblateness suppress orbit flips due to the ELK effect, and generate flips by the ELK-J 2 effect.Moreover, when the stellar oblateness perturbations are of the same order as the octupole perturbations, i.e., κ Ob ∼ ε Oct , the ELK effect is almost entirely suppressed.At nearly the same point, the ELK-J 2 effect reaches its peak for the given ε Oct .However, from a global view, stellar oblateness always reduces the flipping orbits in an aligned configuration.4. For a misaligned system, if the stellar obliquity is small, the system has similar orbit-flipping behaviors as the aligned one.In contrast, if the stellar obliquity is large, the phase space structure changes significantly.
Our results regarding the ELK-J 2 effect are of significance for various astrophysical settings: the migration of planets in the early stage of an exoplanetary system with an oblate stellar and an eccentric exterior companion, which can be a giant planet or a distant star; satellite dynamics around an oblate planet with an eccentric perturber; or secular evolution of the debris disk around a star or a planet accompanied by a distant eccentric massive object.
This work has been supported by the National Natural Science Foundation of China under grant 11872007, and the Fundamental Research Funds for the Central Universities.

Appendix A Octupole-level Equations of Motion
Equations in Newtonian form.The motion of a particle, subject to stellar oblateness and tidal gravitational force from an outer perturber, is governed by The equations are presented in the inertial coordinate system centered at the central star.G is the gravitational constant, r is the position vector of the particle relative to the central star with magnitude, and Ob y and Per y are the perturbing potentials due to the stellar oblateness and the outer perturber's gravity, respectively.We have In the above equations, J 2 and R 0 are the oblateness and radius of the central star, n 0 is the stellar spin axis, and r Per is the position vector of the outer perturber relative to the central star.
Secular equations in vectorial form.We present the secular equations of motion up to octupole order in a vectorial formalism in terms of j and e.We write the vectorial elements j and e in terms of the position r and velocity v: According to first-order perturbation theory, the secular equations of motion can be constructed by calculating the variations due to these averaged perturbing effects separately and then summing them.By substituting Equation (3) into Equations ( 9) and (10), the secular changes in j and e due to the stellar oblateness are given by In the above equations, we have defined a dimensionless time τ = t/t LK , where t LK is the timescale of the LK cycles at the quadrupole order, given by Equation (1).It also takes a form as Note that, although the timescale of the LK-J 2 oscillations deviates from the timescale of the standard LK cycles due to the pericenter advance induced by the stellar oblateness, the scale we choose to deal with does not influence the final results.Finally, taking a sum of them, we obtain the secular equations: e e e e , A 1 2 where j  and e  denote dj/dτ and de/dτ, respectively.
From Equations (A5)-(A11), we find that, in the test-particle limit (n 0 ˆand j Per ˆare fixed), the orbital evolutionary behaviors follow a reflection symmetry.If we replace j with −j, both j  and e  reverse their signs, j j   and e e    -, and one can restore invariance through time reversal, τ → −τ.This operation is equivalent to taking i with π − i, ω with π − ω, and Ω with π + Ω.Moreover, this reflection symmetry is irrelevant to the coplanar configuration between the stellar equator and the perturber's orbit.
Secular equations in Hamiltonian form.Our secular equations in terms of j and e can be easily transformed into Hamiltonian equations based on canonical variables, which are widely used in previous studies of the secular dynamics of the hierarchical triple systems (Kozai 1962;Ford et al. 2000;Lithwick & Naoz 2011;Naoz et al. 2013).For the convenience of confirmation, we convert our equations into Hamiltonian form.
We adopt the normalized canonical Delaunay variables, defined as Note that j = G and j z = H.One can easily derive the relations between the vectorial elements and these canonical variables by expressing j and e in the inertial coordinate system defined in the main body.We have By using the above relations, the double-averaged potentials in Equations (3)-( 5), written in the canonical variables, are ´--+ - In terms of the canonical variables shown in Equation (A13), the dimensionless Hamiltonian is given by Finally, the secular equations are (with respect to the dimensionless time τ  ´-+ ---- ´--

Appendix B Massless Approximation of the Inner Planet
In this paper, the inner planet is assumed to be massless, and thus its effects on the outer orbit and the star's obliquity are ignored.Here, we examine the validity of the two simplifications.
In a hierarchical three-body system, the ratio of the timescales of the secular evolution of the inner and outer orbits depends on the ratio of their angular momenta, i.
The ratio is always much smaller than 1 if the outer perturber is far enough from and more massive than the inner planet.In this case, the outer orbit can be treated as fixed in the inertial space.
Next, we consider the interaction between the angular momenta of the inner orbit L ^^^^T he ratio of the precession timescales of the two unit vectors is determined by L 1 /L 0 .
On the other hand, because the total angular momentum of the inner system L in = L 0 + L 1 is conserved (ignore the outer perturber and the inner planet is treated as a point mass), the two vectors L 0 and L 1 actually precess around L in , as shown in Figure 14.Therefore, alternatively, we can use the angles θ 0 and θ 1 , defined in the picture, to geometrically determine the secular changes in the orientations of the two vectors L 0 and L 1 due to the stellar quadrupole.The angles θ 0 and θ 1 are given by Specifically, we consider such a typical system: a planet moves in an orbit with semimajor axis a = 1 au and eccentricity e = 0.1 around a central star with mass m 0 = m e and radius R 0 = R e .The stellar oblateness J 2 is 1 × 10 −4 , corresponding to 6.49 10 a 5

Figure 1 .
Figure1.Geometry of the hierarchical three-body system.An oblate center star is orbited by a massless inner planet and a massive exterior perturber.The perturber's orbit pole and the stellar spin are assumed to be aligned.

Figure 2 .
Figure2.An example of an orbit flip by the ELK-J 2 effect (left column, red) and comparison to an orbit flip by the ELK effect (right column, green).The blue trajectory in the left column is the numerical integration of the full equations of motion, Equation (A1).In this paper, a fourth-and fifth-order Runge-Kutta method with adaptive step size is used to numerically solve the equations.The two cases of the left and right columns have similar system parameters except the stellar oblateness J 2 , and similar initial conditions except Ω 0 .The system parameters are m 0 = 2 M e , R 0 = 3 R e , J 2 = 10 −3 (J 2 = 0 for the right column), m Per = 50 M J , the particle's mass m = m Per , a = 0.61 au, a Per = 20 au, and e Per = 0.79.With these parameters, the characteristic parameters of the system are κ Ob ; 0.1775 (κ Ob = 0 for the right column) and ε Oct ; 0.0393.The initial conditions of the particle are set as i = 74°, e = 0.87, ω = 0°, and Ω = 0°(180°) for the left (right) column.The total integration time is 1.96 Myr (∼100t LK ).It shows a good agreement between the full equations of motion and the secular model of the octupole-level approximation.

Figure 3 .
Figure3.Phase portraits of a system with κ Ob = 0.1775 at the quadrupolelevel dynamics.From the upper panel, (A) j z = 0.259, (B) j z = 0.136, and (C) j z = 0.105, respectively.The green and red curves are separatrices of the librating (in black) and circulating cycles (in yellow).There are four types of quadrupole-level cycle: librating and circulating LK cycles, and librating and circulating LK-J 2 cycles.

Figure 5 .
Figure 5. Evolution of j z and trajectories in the e i , ( ) space for the triple systems presented in Figure 2. The left column represents the ELK-J 2 effect and the right column represents the ELK case, corresponding to the left and right columns of Figure 2, respectively.The black curves in the lower two planes are contour lines of j e i 1 c o s z 2 = -.

4. 1 .
Phase Space for a Fiducial Example with κ Ob = 0.1, ε Oct = 0.03 First, a specific example with κ Ob = 0.1 and ε Oct = 0.03 is considered.The exploration of its phase space is conducted as follows: the most representative Poincaré sections with a varying energy constant,

Figure 6 .
Figure 6.Contours of the energy constant, 16 15 0 y y y » * /(), in the e i , cos ( ) plane with ω = 0 for a system with κ Ob = 0.1.For clarity, we adopt a different abscissa scale for the e > 0.9 areas.

Figure 7 .
Figure 7. Poincaré surfaces of section at different levels of energy constant presented in the i , cos W ( ) plane for a system specified by κ Ob = 0.1, ε Oct = 0.03.Only the trajectories associated with orbit flips are recorded.The sections in the upper and lower panels are defined by 0, 0  w w = > and 0, 0  w w = < , respectively.In each panel, the green dots characterize the boundaries of the libration islands that allow orbit flips, and i cos D represents the size of the islands.

Figure 8 .Figure 9 .
Figure8.Analytical results: orbit-flipping boundaries in the entire (e 0 , i 0 ) space for a system with κ Ob = 0.1, κ Oct = 0.03 obtained by exploiting the Poincaré surfaces of section.The dominant two regions are designated as I and II, respectively.(Note: although the Poincaré surface of section is based on numerical integrations of the secular equations, we would like to refer to the results acquired in this manner as "analytical" to distinguish them from the numerical exploration in the whole phase space.)

Figure 11 .
Figure11.Minimum flipping inclination due to the ELK effect.Each curve marks the minimum i 0 that induces flips for the given value of ε Oct and different values of κ Ob .Below each curve none of the orbits flip, and above it some orbits do.These curves are obtained by numerically integrating Equations (9) and (10) initialized with e 0 = 0.01, ω 0 = 0, and Ω 0 = π, scanning through the values of i 0 from 0 to π/2 to determine the minimum i 0 that allows flips.We did not scan the initial Ω 0 , because according to the phase space structure at e 0 ∼ 0, the lower limit of flipping i 0 is determined at Ω 0 = π.At this initial e 0 , orbit flips occur in the region dominated by the ELK effect.These curves are symmetric about i 0 = 90°, and only the parts below 90°are depicted.
e., L 1 /L 2 , where -type star, k 0 = 0.1, and a W • is the stellar spin velocity.The relation between a W • and the stellar oblateness J 2 is given by Equation (2).In the stellar quadrupole field,

Figure 15 .
Figure 15.Variations of θ 0 and θ 1 with the mass of the inner planet.