Dynamical Structures under Nonrestricted Hierarchical Planetary Systems with Different Mass Ratios

Secular dynamics have been extensively studied in both the inner and outer restricted hierarchical three-body systems. In the inner restricted problem, the quadrupole-order resonance (i.e., the well-known Kozai resonance) causes large coupled oscillations of eccentricity and inclination when the maximum inclination is higher than 39.2°, and the octupole-order resonance leads to the behavior of orbital flips. In the outer restricted problem, the behavior of orbital flips is due to the quadrupole-order resonance. Secular dynamics under the inner and outer restricted systems are distinctly different. The mass ratio of inner and outer bodies could change the ratio of circular orbital angular momenta β, which significantly influences dynamical structures of the system. But this influence is still unclear. In this paper, we focus on nonrestricted hierarchical planetary systems where β > 1 and investigate the secular dynamics by changing mass ratios. Dynamical structures are systematically explored from four aspects: periodic orbits, secular resonances, orbital flips, and chaos detection. We find that (a) it tends to lead to more bifurcations in the host family of prograde periodic orbits associated with Kozai resonance with smaller β; (b) with the decrease of β, fewer orbits inside the octupole-order resonance can realize flip; (c) for given initial conditions, the forbidden region appears in the retrograde region and becomes larger as β decreases, meaning that the mutual inclination cannot reach a very high value if β is small; and (d) chaotic orbits are distributed in the low-eccentricity, high-inclination region when β > 1.


Introduction
The model of a hierarchical three-body system has been widely applied to diverse astrophysical settings, such as binary asteroids perturbed by the Sun, hot Jupiters with the host star, stars with a black hole, and so on (Naoz 2016).In a hierarchical three-body system, the outer orbit is much farther away from the inner orbit, indicating that the ratio of the semimajor axes of inner and outer orbits is a small quantity.Therefore, the Hamiltonian of the system can be expanded into a power series of the semimajor axis ratio (Harrington 1968(Harrington , 1969)).The Hamiltonian is then double averaged (e.g., von Zeipel transformation) over the inner and outer orbital periods to study the secular dynamics.
When the inner body is a test particle with axisymmetric perturbation from the outer body, the octupole term vanishes, so the Hamiltonian up to the quadrupole order is a good approximation.It is found that the highly inclined test particle would undergo coupled oscillation of eccentricity and inclination at the quadrupole-order approximation (von Zeipel 1910;Kozai 1962;Lidov 1962).This phenomenon is the so-called "von Zeipel-Lidov-Kozai" (ZLK) effect (Lithwick & Naoz 2011;Ito & Ohtsuka 2019).In this situation, secular dynamics are dominated by the quadrupole-order resonance, also called Kozai resonance.The z-component of the inner orbital angular momentum is conserved, so orbital flip cannot happen in this setting (Naoz et al. 2013).
If the inner body is perturbed by an eccentric outer body and masses of the inner binary are not comparable, the octupole order of the Hamiltonian is no longer negligible and needs to be considered.By applying the Hamiltonian up to the octupole order, Naoz et al. (2011) successfully produced hot Jupiters in retrograde orbits.The eccentricity of the perturbed body can be excited to extremely high values, and the inclination can flip from prograde to retrograde and back again because of the long-term modulation of the ZLK cycles.This dynamical behavior is attributed to the eccentric ZLK effect (Lithwick & Naoz 2011;Naoz 2016).In the inner test particle approximation, many studies have been carried out.Katz et al. (2011) derived the conditions for the behaviors of orbital flips and high eccentricity analytically, and Lithwick & Naoz (2011) gave the conditions numerically.By means of the Poincaré sections and the Lyapunov exponents, Li et al. (2014a) found that the eccentric ZLK effect is due to the resonances in the octupoleorder Hamiltonian and chaos occurs with high initial mutual inclination of the inner and outer orbits.In addition, not only the orbits with low eccentricity and high inclination but also the coplanar orbits with high eccentricity can realize flip (Li et al. 2014b).The timescales of the ZLK oscillations at the quadrupole order and of the eccentric ZLK oscillations at the octupole order were given by Antognini (2015).With the help of perturbative treatments, Sidorenko (2018) interpreted the eccentric ZLK effect as a resonance phenomenon by introducing "action-angle" variables.The essence of flipping orbits caused by the eccentric ZLK effect in the test particle approximation is discussed by Lei (2022), showing that flipping orbits are quasi-periodic orbits around the polar periodic orbit from the viewpoint of dynamical system theory and are resonant orbits from the viewpoints of the Poincaré section and perturbation theory.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
When the outer body is a test particle, the behavior of orbital flips occurs mainly by the quadrupole-order resonance (Naoz et al. 2017).The critical argument for resonance-caused orbits to flip is the longitude of the ascending node of the outer orbit, centered at ±90° (Naoz et al. 2017;Zanardi et al. 2017;Vinson & Chiang 2018).This indicates that the essences of orbital flips are different between the inner and outer restricted problems, although the behavior of orbital flips in the outer restricted system is similar to that in the eccentric ZLK effect.Furthermore, when the Hamiltonian is up to the hexadecapole order, the inverse Kozai resonance occurs (Vinson & Chiang 2018).
For the nonrestricted hierarchical three-body systems, where the masses of the inner and outer bodies cannot be neglected, dynamics are different from the restricted case in some way.The z-component of the inner orbital angular momentum is no longer conserved, and even at the quadrupole-order approximation the behavior of orbital flips may occur (Naoz et al. 2013).Secular dynamics under the parameters, including semimajor axes, mutual inclination, masses, and eccentricity of the outer planet, are revealed numerically by Teyssandier et al. (2013), showing that dynamics (extreme eccentricity and orbital flips) beyond the test particle limit are different from those in the restricted systems.The flipping orbits caused by the octupole-order resonance are around resonant orbits without flips, due to the asymmetry of the Hamiltonian in the nonrestricted systems (Huang & Lei 2022).Hansen & Naoz (2020) surveyed the stationary points at both the quadrupoleorder approximation and the octupole-order approximation to find possible secular dynamics in the nonrestricted planetary systems.They found that the mass ratios have an influence on the apsidal alignment (Hansen & Naoz 2020).Hamers (2021) presented a semianalytical approach for properties of the ZLK oscillations at the quadrupole-order approximation.The quadrupole-order and octupole-order resonances in the nonrestricted hierarchical systems, where the outer body is more massive than the inner one, are studied by Lei & Huang (2022) from the viewpoint of perturbative theory.The behavior of the eccentric ZLK effect combined with tidal dissipation has been applied to the formation of hot Jupiters (Naoz et al. 2011;Petrovich 2015;Dawson & Johnson 2018).More applications of the eccentric ZLK effect can be found in Naoz (2016) and Shevchenko (2016).
However, the question of how the dynamical structures are affected by the masses of the inner and outer planets remains unanswered, since the dynamics are very different between the inner restricted system and the outer restricted one.The mass ratio could change the ratio of angular momenta of the inner and outer circular orbits, which influences the dynamics of the system.In this paper, we carry out a systematic study about the influence of planetary masses on secular dynamics in nonrestricted planetary systems.We discuss the cases with the ratio of circular orbital angular momenta greater than 1.First, the Hamiltonian in the hierarchical system up to the octupole order in the semimajor axis ratio is formulated.Second, characteristic curves of periodic orbits in the quadrupole-order and octupole-order resonances are presented in the parameter spaces to have an overview of dynamics under different mass ratios.Third, distributions of varieties of resonances are mapped out under different mass ratios.Fourth, flipping regions are explored.By comparing resonant regions and the flipping regions, the regions of orbital flips caused by resonances can be well recognized.Finally, a normalized second-derivative-based index is introduced for chaos and separatrix detection.Dynamics can be well understood thanks to dynamical maps, which have excellent correspondence to the distributions of periodic orbits and resonant regions.
The structure of this paper is organized as follows: In Section 2, the Hamiltonian model is presented.Periodic orbits are shown in Section 3. Resonant regions and flipping regions are explored in Sections 4 and 5, respectively.In Section 6, dynamics affected by different mass ratios are revealed according to the chaos and separatrix detection.Our conclusions are summarized in Section 7.

Hamiltonian Model
In a nonrestricted hierarchical planetary system, the outer planet with mass m 2 moves around the inner binary with masses m 0 and m 1 on a much wider orbit (m 1 , m 2 = m 0 ).The dynamical model is usually formulated in the Jacobi coordinate frame, where the z-axis is along the direction of the total angular momentum G tot and the x-y plane is coincident with the invariant plane of the system (Naoz et al. 2013).Under the coordinate system, the semimajor axis, eccentricity, argument of pericenter, longitude of ascending node, and orbital inclination are denoted by a, e, ω, Ω, and i, respectively.For convenience, subscripts "1" and "2" represent the inner and outer orbits, respectively.The mutual inclination between the inner and outer orbits is denoted by i tot .
Due to the hierarchical configuration, the Hamiltonian is usually truncated at octupole order and doubly averaged over orbital periods of the inner and outer binaries to get secular dynamical models (Harrington 1968(Harrington , 1969)).The resulting normalized long-term Hamiltonian is given by (Krymolowski & Mazeh 1999;Ford et al. 2000;Naoz et al. 2013 To derive equations of motion, Delaunay's variables are introduced as where β is Here G 1 and G 2 are normalized orbital angular momenta of the inner and outer orbits.Given that in a planetary system the masses of planets are negligible compared to the stellar mass, Equation so for a given semimajor axis ratio, mass ratio influences secular dynamics by affecting the ratio of circular orbital angular momenta β.The relation about longitudes of the ascending node Ω 1 − Ω 2 = π always holds (Naoz et al. 2013), so longitudes of the ascending node are absent in the Hamiltonian.However, this does not indicate that H 1 and H 2 are constants (Naoz et al. 2013).
The equations of motion can be derived from the Hamiltonian canonical relations, and the mutual inclination i tot can be derived from the relation of angular momenta (Naoz et al. 2013), The dynamical system has two degrees of freedom, depending on the total angular momentum G tot .From Equations (1)-( 3) and (6), we find that the normalized Hamiltonian  is independent of m 2 and m 1 for planetary systems, and the planetary mass ratio m 2 /m 1 only contributes to G tot .To understand how dynamics are affected by the masses of planets, we could take different planetary mass ratios.In this paper we discuss the dynamics transition in the cases of β > 1.We set m 0 = 1M e and m 1 = 1M J , where M e represents the mass of the Sun and M J represents the mass of the Jupiter.The masses of m 2 are set as 40M J , 20M J , 10M J , 5M J , 2.5M J , 1.1M J , 0.9M J , and 0.3M J , corresponding to β = 168, 84.8, 42.6, 21.3, 10.7, 4.70, 3.85, and 1.28.Other parameters of systems are taken from the exoplanetary system HD 4113, which are a 1 = 1.298 au and a 2 = 23.72618au (Tamuz et al. 2008;Cheetham et al. 2018).In our previous work (Huang & Lei 2022), we have shown numerically that, with different values of G tot , dynamical structures are qualitatively similar for nonrestricted planetary systems.Thus, in this work, the total angular momentum is taken as which corresponds to e 1 = e 2 ≡ 0.3 and i tot ≡ 90°in Equation (6).The parameters of systems can be found in Table 1.Furthermore, when the initial values of G 1 and i tot are given, the value of G 2 can be solved by ), there is no physical solution ( i.e., forbidden regions).Secular dynamics in a hierarchical planetary system are dominated by the quadrupole-order resonance and the octupole-order resonance (Lei & Huang 2022).According to Equation (1), the quadrupole-order resonance (i.e., Kozai resonance) is centered at g 1 = ± π/2 (Kozai 1962).When Kozai resonance happens, g 1 is librating around ±π/2, and g 2 is circulating.The octupole-order resonance is found to be of g 1 + g 2 libration or g 1 − g 2 libration, depending on the orbital inclination, while 3g 1 ± g 2 resonance cannot be found in our numerical exploration.For both prograde and retrograde cases, the critical argument of the octupole-order resonance is generally denoted by s 1 1 , librating around 0 or π.Here sign(x) is equal to 1 when x > 0 and equal to −1 when x < 0. Based on these analyses, when exploring dynamics dominated by the quadrupole-order term in the Hamiltonian, we take the initial values of g 1 and g 2 fixed at g 1 = π/2 and g 2 = π/2 for initial antialigned configurations.In the octupole-order exploration, the initial value of g 1 is set at g 1 = 0, which is a saddle point of Kozai resonance, and the initial values of g 2 are taken as g 2 = 0 for the case σ = 0 and g 2 = π for the case σ = π.
In the following sections, we numerically explore the quadrupole-order resonance and octupole-order resonance from four aspects: periodic orbits, secular resonances, orbital flips, and chaos detection.We use Runge-Kutta-Fehlberg (RKF78) to integrate equations of motion represented by Equation (5) (Fehlberg 1969).The tolerance adopted for integrations is 10 −12 .Dynamical structures under different mass ratios are discussed.

Periodic Orbits
In this section, dynamical structures under different mass ratios are discussed from the aspect of periodic orbits.Let Φ(X 0 , t) be the phase flow of an orbit at time t with initial state From the viewpoint of the Poincaré section, periodic orbit corresponds to the center of the libration island.Thus, we take the center arising in the Poincaré section as an approximate guess of the periodic orbit, denoted as ¢ X 0 , and the associated orbital period is denoted by ¢ T .The initial guesses are further corrected by Substituting Equation (10) into Equation (9) and expanding Equation (9) at ¢ ¢ ( ) X T , 0 to the first order, we can get an iterative equation for determining the periodic orbit (Deprit & Henrard 1967): Here I is the identity matrix.The partial derivatives can be derived by the central difference.Observing from Equation (11), there are five parameters to be solved, but there are only four equations.However, a periodic orbit is a resonant center in the Poincaré section, and the initial angles can be fixed at the resonant center, meaning that Δg 1 = 0 and Δg 2 = 0.In this manner, the periodic orbit can be determined by iterating Equation (11).
To produce a family of periodic orbits, the predictioncorrection algorithm is introduced after getting an initial periodic orbit.Let X 0 1 and T 1 be initial states and the period of the known periodic orbit.We suppose that there is a nearby periodic orbit ( ) X T , 0 2 2 , with the relation The prediction equation is derived by substituting Equation (12) into Equation (9) and expanding Equation (9) at ( ) X T , 0 1 1 up to the first order (Deprit & Henrard 1967): Please note that Δg 1 = 0 and Δg 2 = 0 are fixed in Equation (13).When a small quantity ΔG 2 is given, other small quantities, ΔG 1 and ΔT, are solved from Equation (13).
Because Equation ( 13) is a linear approximation, the solutions ΔX and ΔT are approximations, denoted as D ¢ X and D ¢ T .The corresponding approximate values are taken into Equation (11) to iterate in order to achieve an accurate periodic orbit.By giving a small change ΔG 2 and applying the prediction-correction algorithm, a family of periodic orbits can be generated.
The stability of the periodic orbit is determined with the help of the monodromy matrix M, where M is the state transition matrix taken at one period of the periodic orbit.For our autonomous Hamiltonian system of two degrees of freedom, there are two reciprocal pairs of the eigenvalues for M, and one pair of them has a value of 1.If the other pair of the eigenvalues λ and λ −1 satisfy |λ + λ −1 | < 2, the periodic orbit is stable (Lara et al. 2007).

Kozai Resonance
The stable periodic orbits corresponding to Kozai resonance under different mass ratios are shown in the i tot -e 1 plane in Figure 1.Stable periodic orbits in systems with β = 168, 84.8, 42.6, and 21.3 are denoted by black, red, blue, and green dots, respectively, in the top panel, and those in systems with β = 10.7,4.70, 3.85, and 1.28 are denoted by black, red, blue, and green dots, respectively, in the bottom panel.From Figure 1, we can see that (a) periodic orbits corresponding to Kozai resonance happen within the inclination interval 39°-141°, which is consistent with Kozai (1962); (b) the family of Kozai prograde periodic orbits has bifurcation compared to retrograde ones, and it tends to have more bifurcation in the host family of prograde periodic orbits with lower β; and (c) when β 10.7, retrograde periodic orbits only distribute in high-eccentricity regions and gradually disappear as β decreases.
We take β = 168 as an example to further explore bifurcations in prograde configurations, as shown in Figure 2. In the top left panel, periodic orbits of the full octupole-order Hamiltonian correspond to the characteristic curve with different colors, where different colors represent different levels of the Hamiltonian, and stationary points of the quadrupole-order Hamiltonian defined in Equation (4) in Hansen & Naoz (2020) are denoted by a black curve.By comparing the two curves, it is known that bifurcations are induced by the octupole-order term.Furthermore, Poincaré sections at different levels of the Hamiltonian are shown in the remaining three panels.The libration centers are marked in Poincaré sections, and the corresponding periodic orbits are marked on the characteristic curve.According to Poincaré sections, the occurrence of bifurcations is due to the split of Kozai libration islands.When the Hamiltonian is at =  5, there is only one Kozai resonant island.As the Hamiltonian is increased to =  7.6, the Kozai resonant island splits into two islands, corresponding to bifurcations in the prograde periodic orbits.When the Hamiltonian is further increased up to =  17, one island disappears and there remains one Kozai resonant island.

The Octupole-order Resonance
Stable and unstable periodic orbits corresponding to the octupole-order resonance centered at σ = 0 are shown in the i tot -e 1 plane in Figure 3. Unstable periodic orbits are denoted by red dots.In the top panel, stable periodic orbits in systems with β = 168, 84.8, 42.6, and 21.3 are denoted by black, magenta, blue, and green dots, respectively.In the bottom panel, periodic orbits in systems with β = 10.7,4.70, 3.85, 1.28 are denoted by black, magenta, blue, and green dots, respectively.From the viewpoint of Poincare sections, stable periodic orbits correspond to resonant centers and unstable periodic orbits correspond to saddle points.We can observe from Figure 3 that (a) generally speaking, there are two families of periodic orbits associated with the octupole-order resonance centered at σ = 0, distributed in the low-and high-eccentricity regions; (b) the families of periodic orbits are more symmetric about i tot = 90°with a larger β; and (c) as β decreases, periodic orbits gradually disappear in the retrograde region.
For the octupole-order resonance centered at σ = π, stable and unstable periodic orbits are shown in the i tot -e 1 plane in Figure 4.It is observed that when β > 10.7, structures of the family of periodic orbits where the stable ones are distributed in the intermediate-eccentricity and near-polar regions are similar to that in octupole-order resonance centered at σ = 0, but stable and unstable regions are reversed.This is because if a resonant island is centered at 0, the corresponding saddle point is at π, and vice versa.Besides, two families of periodic orbits are distributed in the low-eccentricity, low-inclination region.With a larger β, distribution of periodic orbits is more symmetric about i tot = 90°.As β decreases, retrograde periodic orbits gradually disappear.Moreover, although it seems that distributions of periodic orbits in the bottom panel differ a lot from .Characteristic curves of stable and unstable periodic orbits corresponding to the octupole-order resonance centered at σ = 0 under different mass ratios.The initial angles are set at g 1 = 0, g 2 = 0.The stable orbits are denoted by the colors shown in the legends, and the unstable periodic ones are all denoted by red dots.
those in the top panel, we will show in Section 4 that resonant islands associated with these periodic orbits in the bottom panel are small.
In the following section, resonant regions are discussed to further reveal dynamics under different mass ratios.

Secular Resonances
In this section, distributions of secular resonances under different mass ratios are given numerically.When integrating the equations of motion represented by Equation (5), if the critical argument is always librating around a certain resonant center, the initial condition is recorded.Kozai resonance and the octupole-order resonances are explored in this section.

Kozai Resonance
If the critical argument g 1 is always librating around g 1 = π/2, the orbit is inside Kozai resonance.Nonresonant orbits are denoted by pure blue dots in the i tot -e 1 plane in Figure 5, and resonant orbits are denoted by different colors that correspond to different values of the Hamiltonian.Forbidden regions are shown in white, where there is no physical solution for the given total angular momentum.The black dots stand for initial conditions of stable Kozai periodic orbits (i.e., Kozai resonant centers).
In general, ignoring forbidden regions, dynamical distribution of Kozai resonance is similar under different mass ratios.Kozai resonance is distributed within mutual inclination between 39°and 141°, in agreement with the previous classic works (Kozai 1962;Naoz 2016;Shevchenko 2016).The excursion of planets can be observed according to the same level of the Hamiltonian.Although contour lines of the Hamiltonian could cross from prograde to retrograde, we will show in Section 6 that, due to the separatrix at i tot = 90°, the behavior of orbital flips can hardly happen.
When β = 21.3, the forbidden region arises in the retrograde region.As β decreases, the forbidden region grows, resulting in the absence of retrograde periodic orbits and the reduction of the Kozai resonant area.

The Octupole-order Resonance
When the octupole-order resonance happens, the critical argument s is librating around σ = 0 or σ = π.The octupole-order resonance centered at σ = 0 is shown in Figure 6, and the resonance centered at σ = π is shown in Figure 7.In Figures 6 and 7, nonresonant orbits are denoted by pure blue, and resonant regions are denoted by different colors where the color bars represent values of the Hamiltonian.The white region is the forbidden zone.The corresponding stable periodic orbits (i.e., resonant centers) are denoted by black dots, and unstable periodic orbits (i.e., saddle points) are denoted by red dots.As expected, resonant centers are surrounded by resonant orbits, and saddle points are distributed outside the resonant regions.With a larger β, which means that the configuration is closer to the inner restricted problem, dynamical structure is more symmetric about i tot = 90°.
We first discuss secular resonance centered at σ = 0.In general, this resonance is distributed in the low-and higheccentricity regions according to Figure 6.Furthermore, by observing the same level of the Hamiltonian, it is found that if the initial e 1 is greater than a critical value e 1,c , a near-coplanar orbit can be inside resonance for both prograde and retrograde configurations when β 21.3.As β decreased from 21.3 to 1.28, the resonant regions are gradually forced into the prograde zone owing to the growth of the forbidden region, and the high-eccentricity resonant region gradually moves to the low-eccentricity area.
Then, the octupole-order resonance centered at σ = π is explored, as shown in Figure 7.This resonance is distributed in the intermediate-eccentricity, high-inclination region and in the low-eccentricity, low-inclination region.With the decrease of β, retrograde resonant regions gradually disappear because of the occurrence of the forbidden region, and the resonant area becomes smaller as well.
Secular resonance may cause orbital flips when it is centered near i tot = 90°for the nonrestricted three-body system, where the outer planet is more massive than the inner one (Huang & Lei 2022).Especially for the inner restricted problem, orbital flips can be interpreted as a resonance phenomenon (Sidorenko 2018), and this is due to the octupole-order resonance centered right at i tot = 90°(Lei 2022).However, for the outer restricted problem, orbital flips are due to the quadrupole-order resonance centered at i tot = 90° (Naoz et al. 2017).In the next section, dynamics of systems with β > 1 are studied from the viewpoint of orbital flips.

Orbital Flips
Orbital flips happen if the mutual inclination i tot can switch between prograde and retrograde along the evolution (Naoz 2016).There are two main kinds of flipping orbits, including regular and chaotic flipping orbits.The regular flipping orbits can be further divided into resonant flipping orbits and circulate flipping orbits.When the outer planet is more massive than the inner one, the majority of flipping orbits are resonant orbits (Huang & Lei 2022).In this section, orbital flips are explored under different mass ratios to reveal the influence of planetary mass on dynamics.By numerical integration over more than 2000 Kozai cycles, orbital flips associated with Kozai resonance and the octupole-order resonance are studied respectively. .Characteristic curves of stable and unstable periodic orbits corresponding to the octupole-order resonance centered at σ = π under different mass ratios.The initial angles are set at g 1 = 0, g 2 = π.The stable ones are denoted by the colors shown in the legends, and the unstable periodic ones are all denoted by red dots.

Orbital Flips Associated with Kozai Resonance
In the i tot -e 1 plane, initial values that can generate flipping orbits are denoted as green dots, as shown in Figure 8.For those that cannot generate flipping orbits, the initial values of orbits inside Kozai resonance are denoted as red dots, and those outside Kozai resonance are denoted as blue dots.Stable periodic orbits are marked by black dots.The forbidden region is shown in white.Here the range of i tot is set between 39°and 141°, where Kozai resonance happens.
In general, flip regions only occupy a small portion of the feasible region, and the distribution of flipping orbits under different mass ratios is similar.Only a few Kozai resonant orbits distributed around i tot = 90°can realize flip.When β 10.7, a few orbits outside resonance can realize flip as well.In Section 6 we will show that the flipping orbits outside resonance are chaotic.

Orbital Flips Associated with the Octupole-order Resonance
The behavior of orbital flips caused by the eccentric ZKL effect is referred to as a resonance phenomenon in the nonrestricted planetary systems (Huang & Lei 2022).The critical argument of the resonance is s   degenerated to the inner restricted problem, where the behavior of orbital flips is due to the octupole-order resonances centered exactly at i tot = 90° (Lei 2022).Comparing the flipping regions to the resonant regions and periodic orbits in Figures 9 and 10, it is found that with a larger β (e.g., β = 168), the problem is more similar to the inner restricted case, structures are more symmetric, and the near-polar periodic orbits can be closer to i tot = 90°.Therefore, the majority of resonant orbits can realize flip when β is large.Some near-polar orbits that are outside resonance can realize flip as well.As β decreases, the nearpolar periodic orbits gradually deviate from the location of i tot = 90°, so inside the libration islands more and more resonant orbits that are around the resonant centers cannot realize flip.When β is decreased to or less than 4.70, the majority of flipping orbits are distributed near the polar region and the majority of resonant orbits cannot flip.We will further  show in Section 6 that the near-polar flipping orbits that are outside resonance are of circulation and flipping orbits in the low-eccentricity, high-inclination area are chaotic.
In conclusion, when the outer planet is more massive than the inner one, the majority of orbital flips are due to the octupole-order resonance.Flipping orbits inside resonance centered at σ = 0 are distributed in the low-and higheccentricity regions, and those inside resonance centered at σ = π are distributed in the intermediate-eccentricity regions.As β decreases, less resonant orbits can realize flip and the circulating flipping orbits that are distributed near the polar region take a larger portion.

Chaos Detection
In this section, chaos detection is applied to the parameter space of initial conditions under different mass ratios.To detect resonant regions, separatrices, and chaotic seas, the index D   D , which is based on the second derivative of orbital amplitudes, is employed (Daquin & Charalambous 2023).The normalized indicator D   D is given by where D(G 1 ; t) and D(G 2 ; t) are the orbital diameters, measured in the Poincaré sections  with the given initial states X 0 and integration time t, When exploring the quadrupole-order resonance, D(G 1 ; t) and D(G 2 ; t) are measured in the Poincaré section defined by g 2 = π/2, and defined by g 1 = 0 when exploring the octupoleorder resonance.Equation ( 14) can be calculated by central differences (Daquin & Charalambous 2023).The values of D   D are at high levels for chaotic orbits, separatrices, and periodic orbits and at low levels for circulate and quasi-periodic orbits (Daquin & Charalambous 2023).Thus, dynamical maps represented by D   D are able to detect resonant centers, secondary resonances, and chaos.We also employ the fast Lyapunov indicator (FLI) to recognize chaos, which can be found in the Appendix.

Dynamical Maps Associated with Kozai Resonance
The dynamical maps associated with Kozai resonance under different mass ratios are shown in Figure 11, and the periodic orbits are marked by black dots.As can be seen, the resonant centers, separatrices, and chaotic regions can be well recognized from the D   D maps.It is observed from Figure 11 that (a) the resonant centers highlighted in bright yellow in the D   D maps are coincident with the periodic orbits; (b) the separatrices distributed in i tot = 39°, i tot = 90°, and i tot = 141°show the boundaries of Kozai resonance; (c) there is another separatrix inside the prograde resonant region, separating two families of bifurcated periodic orbits; and (d) the low-eccentricity regions highlighted in yellow are filled with chaotic orbits.

Dynamical Maps Associated with Octupole-order Resonance
Dynamical maps associated with the octupole-order resonance centered at σ = 0 and σ = π are shown in Figures 12 and  13, respectively.The corresponding stable periodic orbits are marked by black dots, which are represented by bright yellow dots in the D   D maps as well.Furthermore, the resonant regions surrounded by separatrices can be easily identified according to the bright yellow lines.
It is found from Figures 12 and 13 that (a) compared to the resonant distribution in Section 4.2, the octupole-order resonant

Conclusions
In this work, dynamical structures are explored numerically under different mass ratios in the nonrestricted planetary systems with β > 1.The Hamiltonian up to the octupole order in semimajor axis ratio is adopted and then averaged over the periods of the inner and outer orbits to explore the secular effects.The variations of dynamics affected by mass ratios are revealed by the periodic orbits, resonant regions, flipping regions, and dynamical maps, which are in good agreement with one another.The features can be summarized as follows.With a larger β, the problem is closer to the inner restricted case.The quadrupole-order resonance (i.e., ω 1 resonance) is distributed in the region of high mutual inclination (39°-141°).There is a bifurcation in the family of prograde periodic orbits associated with the quadrupole-order resonance, and the family of retrograde periodic orbits associated with the quadrupoleorder resonance has no bifurcations.The appearance of bifurcation is due to the secondary resonance inside Kozai resonance.
The octupole-order resonance (i.e., -( ) ( ) i g i g sign cos sign cos 2 2 1 1 resonance) centered at σ = 0 is distributed in the low-and high-eccentricity regions.The octupole-order resonance centered at σ = π is distributed in the intermediate-eccentricity region and in the low-eccentricity, low-inclination regions.The behavior of orbital flips is mainly due to the octupole-order resonance centered near i tot = 90°.
As β decreases, dynamical structures gradually change.First, it tends to lead to more bifurcations in the primary family of prograde periodic orbits associated with the quadrupole-order resonance.Second, the location of the octupole-order resonance gradually deviates from i tot = 90°, so more and more resonant orbits around the resonant centers cannot realize flip.The portion of circulating flipping orbits is increased, distributed around i tot = 90°.Moreover, with a smaller β, the forbidden region is larger, meaning that the mutual inclination of two orbits cannot be very high for the given initial conditions.Because of the forbidden region, periodic orbits, resonant orbits, and flipping orbits gradually disappear in the retrograde region.
The behavior of chaos occurs when the eccentricity of the inner orbit is low and the mutual inclination is high for β > 1.

Figure 1 .
Figure 1.Characteristic curves of stable periodic orbits associated with Kozai resonance under different mass ratios in the i tot -e 1 plane.The initial angles are fixed at g 1 = π/2 and g 2 = π/2.

Figure 2 .
Figure2.Top left panel: periodic orbits associated with Kozai resonance in the system of β = 168 of the full octupole Hamiltonian (colored curve), and fixed points defined in Equation (4) inHansen & Naoz (2020) under the quadrupole-order Hamiltonian (black curve, for which e 2 is set at 0.3).The last three panels: Poincaré sections defined by g 2 = 90°at different levels of the Hamiltonian.The resonant centers are marked by red asterisks with letters from "A" to "D."

Figure 4
Figure 4. Characteristic curves of stable and unstable periodic orbits corresponding to the octupole-order resonance centered at σ = π under different mass ratios.The initial angles are set at g 1 = 0, g 2 = π.The stable ones are denoted by the colors shown in the legends, and the unstable periodic ones are all denoted by red dots.

Figure 5 .
Figure 5. Dynamical structures associated with Kozai resonance under different mass ratios.The regions denoted by different colors stand for orbits inside Kozai resonance, and the color bars represent different values of the Hamiltonian.The pure blue regions stand for orbits outside Kozai resonance.The black dots represent stable periodic orbits associated with Kozai resonance.Forbidden regions are shown in white.The initial angles are fixed at g 1 = π/2 and g 2 = π/2.

Figure 6 .
Figure 6.Dynamical structures associated with the octupole-order resonance centered at σ = 0 under different mass ratios.The regions denoted by different colors stand for orbits inside resonance, and the color bars represent different levels of the Hamiltonian.The pure blue regions stand for orbits outside resonance.The black dots represent stable periodic orbits, and the red dots represent unstable periodic orbits.The white area stands for the forbidden region.The initial angles are fixed at g 1 = 0 and g 2 = 0.
, librating around σ = 0 or σ = π.The results of orbital flips are shown in Figures 9 and 10, respectively.The initial values of flipping orbits are denoted as green dots.The resonant orbits without flipping are denoted as red dots, and the orbits outside resonance without flipping are denoted as blue dots.Stable and unstable periodic orbits are marked by black and yellow dots, respectively.Initial angles are set as g 1 = 0 and g 2 = 0 in Figure 9 and as g 1 = 0 and g 2 = π in Figure 10.If β → ∞ , the problem

Figure 7 .
Figure 7. Same as Figure 6, but the initial angles are fixed at g 1 = 0 and g 2 = π, representing dynamical structures associated with the octupole-order resonance centered at σ = π.

Figure 8 .
Figure 8. Dynamical structures at initial values of g 1 = π/2 and g 2 = π/2 under different mass ratios.The green dots represent flipping orbits.The red dots are Kozai resonant orbits that cannot flip.The blue dots represent orbits that are outside Kozai resonance and cannot flip.Stable periodic orbits are represented by black dots.The white area stands for the forbidden region.

Figure 9 .
Figure 9. Dynamical structures at initial values of g 1 = 0 and g 2 = 0 under different mass ratios.The green dots represent flipping orbits.The red dots are the octupole-order resonant orbits that librate around σ = 0 and cannot flip.The blue dots represent orbits that are outside resonance and cannot flip.Stable periodic orbits are represented by black dots, and unstable ones are denoted by yellow dots.The white area stands for the forbidden region.

Figure 10 .
Figure 10.Same as Figure 9, but the initial angles are set at g 1 = 0 and g 2 = π, representing orbital flips about octupole-order resonance centered at σ = π.

Figure 11 .
Figure 11.Dynamical maps of the normalized second-derivative-based index D   D corresponding to the quadrupole-order resonance.The initial angles are fixed at g 1 = π/2 and g 2 = π/2.Black dots correspond to the stable periodic orbits associated with Kozai resonance.

Figure 12 .
Figure 12.Dynamical maps of the normalized second-derivative-based index D   D associated with the octupole-order resonance centered at σ = 0.The initial angles are fixed at g 1 = 0 and g 2 = 0. Black dots stand for the stable periodic orbits associated with the octupole-order resonance centered at σ = 0.

Figure 13 .
Figure 13.Dynamical maps of the normalized second-derivative-based index D   D associated with the octupole-order resonance centered at σ = π.The initial angles are fixed at g 1 = 0 and g 2 = π.Black dots stand for the stable periodic orbits associated with the octupole-order resonance centered at σ = π.

Figure 14 .
Figure14.The FLI maps under different mass ratios with the initial angles of g 1 = π/2 and g 2 = π/2.The chaotic orbits are marked by yellow dots.The white area stands for the forbidden region.

Figure 15 .
Figure15.The FLI maps under different mass ratios with the initial angles at g 1 = 0 and g 2 = 0.The chaotic orbits are marked by yellow dots, and the regular orbits are marked by blue dots.The white area stands for the forbidden region.

Table 1
System Parameters under the Configurations of m 0 = 1M e , m 1 = 1M J , and a 1 /a 2 = 0.0547