Dynamical Structures Associated with High-Order and Secondary Resonances in the Spin–Orbit Problem

In our solar system, spin–orbit resonances are common under Sun–planet, planet–satellite, and binary asteroid configurations. In this work, high-order and secondary spin–orbit resonances are investigated by taking numerical and analytical approaches. Poincaré sections as well as two types of dynamical maps are produced, showing that there are complicated structures in the phase space. To understand numerical structures, we adopt the theory of perturbative treatments to formulate resonant Hamiltonian for describing spin–orbit resonances. The results show that there is an excellent agreement between analytical and numerical structures. It is concluded that the main V-shape structure arising in the parameter space (θ̇,α) is sculpted by the synchronous primary resonance, those minute structures inside the V-shape region are dominated by secondary resonances and those structures outside the V-shape region are governed by high-order resonances. Finally, the analytical approach is applied to binary asteroid systems (65803) Didymos and (4383) Suruga to reveal their phase-space structures.


INTRODUCTION
Spin-orbit resonance is a common phenomenon in the Universe and it takes place if the rotational period of one body and its orbital period around the central object are commensurable (Goldreich & Peale 1966;Peale 1977).For example, the planet Mercury is located inside the 3:2 spin-orbit resonance (Lemaitre et al. 2006;Peale & Gold 1965), and our Moon is locked inside a synchronous (1:1) spin-orbit resonance, resulting in the fact that its one side is permanently facing our Earth (Peale 1969).Additionally, chaotic rotations are observed for minor objects in our Solar system, such as Saturn's satellite Hyperion (Wisdom et al. 1984), Nyx and Hydra (Correia et al. 2015;Showalter & Hamilton 2015).In recent years, it is found that spin-orbit resonance commonly exists in binary asteroid systems (Pravec et al. 2016;Naidu & Margot 2015;Scheeres et al. 2006;Pravec et al. 2019).
The conventional model of spin-orbit resonance is formulated under the assumption that the two bodies move Corresponding author: Hanlun Lei leihl@nju.edu.cnaround their barycenter in fixed Keplerian orbits (Goldreich & Peale 1966;Celletti 1990a,b).Under the classical model, spin-orbit resonance has been extensively explored.For Saturn's satellite Hyperion, Wisdom et al. (1984) studied the mechanism for the onset of chaos through the overlap of first-order resonances (Chirikov 1979) and they predicted that there is a large-scale chaotic behavior if the asphericity parameter α satisfies Under a nearly-integrable Hamiltonian model, Celletti & Chierchia (2000) numerically explored the stability of periodic orbits related to spin-orbit resonances by checking the eigenvalues of the monodromy matrix.Flynn & Saha (2005) formulated second-order Hamiltonian model for describing spin-orbit primary resonances (not secondary resonances) by taking Lie-series transformation theory.Jafari Nadoushan & Assadian (2015, 2016) provided explicit expressions for primary spin-orbit resonance width in terms of asphericity and eccentricity and updated the Chirikov criterion for the onset of chaos.
In addition, they thoroughly explored the effects of asphericity and eccentricity upon the dynamical struc- If the orbit and rotation of a binary are coupled, it becomes the spin-orbit coupling problem.In this respect, Naidu & Margot (2015) simulated the coupled spin and orbital motions for binary asteroid systems using Poincaré surfaces of section.Hou et al. (2017) provided explicit expressions to compute the mutual potential, force and torque between two rigid bodies with arbitrary shapes and mass distributions.Under the ellipsoid-ellipsoid model, Hou & Xin (2017) presented modified Hamiltonian models for describing spin-orbit, spin-spin and spin-orbit-spin resonances, showing that the spin-orbit resonance center may change with the mass ratio and mutual distance and the center of spinspin and spin-orbit-spin resonances may change with the rotation states of the binary asteroid.Recently, Jafari-Nadoushan (2023) offered a criterion for dynamical closeness of asteroid pair based on the difference between rotational and orbital angular momenta, and the author explored the dependence of phase-space structures upon the semimajor axis, eccentricity, mass ratio of the secondary to the primary and equatorial elongation of the secondary.
Almost all the works mentioned above concentrate on primary spin-orbit resonances.However, if the frequency of small-amplitude oscillations about the primary resonance and the orbital frequency are commensurable, secondary spin-orbit resonances may happen.The existence of secondary resonances make the phasespace structures be more complicated.In this respect, Wisdom (2004) developed a perturbative model to study the secondary 3:1 spin-orbit resonance of Enceladus and discussed the rate of tidal dissipation for Enceladus inside this secondary resonance.Gkolias et al. (2016Gkolias et al. ( , 2019) ) took advantage of Lie-series transformation theory to formulate analytical approximations for those loworder secondary resonances (1:1, 2:1 and 3:1) bifurcating from the synchronous primary resonance.For the bifurcation diagram of secondary 1:1 resonance, the analytical (normal form) solution can agree with numerical result in the range of α ∈ [0, 1.2], but starts to diverge when α is greater than 1.2 (Gkolias et al. 2019).
It is not difficult to see that the Hamiltonian model based on Lie-series transformation adopted by Flynn & Saha (2005) can describe high-order primary resonances outside the primary island but cannot describe secondary resonances inside the primary island, while the Hamiltonian model based on Lie-series method adopted by Gkolias et al. (2016Gkolias et al. ( , 2019) ) can be used to study sec-ondary resonances but cannot be used to study those high-order spin-orbit resonances.Thus, it is required a unified resonant Hamiltonian model for describing highorder and secondary resonances in the spin-orbit problem to explain phase-space structures.
In this work, two types of dynamical maps (the fast Lyapunov indicator and the normalized secondderivative-based index) are produced in the ( θ, α) space.There are distinguishable structures, governed by synchronous primary resonance, high-order resonances outside the primary resonance and secondary resonances inside the primary resonance.To understand these numerical structures, the theory of perturbative treatments (Henrard & Lemaitre 1986;Henrard 1990) is taken to formulate a unified resonant Hamiltonian model, which can be used to study both the high-order and secondary resonances analytically.Results show that there is an excellent correspondence between analytical structures arising in phase portraits and numerical structures arising in dynamical maps.
The remaining part of this work is organized as follows.In Section 2, the Hamiltonian function of the spinorbit problem is briefly introduced.In Section 3, two types of dynamical maps are produced and Section 4 takes advantage of the theory of perturbative treatments to formulate the resonant Hamiltonian model.Applications to binary asteroid systems are presented in Section 5 and conclusions are summarized in Section 6.

HAMILTONIAN MODEL
The spin-orbit problem states that a triaxial satellite moves around the central planet on a fixed Keplerian orbit with semimajor axis a and eccentricity e (Murray & Dermott 1999).To simplify the problem, we consider a planar configuration where the satellite's spinaxis (along its shortest physical axis) is assumed to be normal to its orbital plane.The principal moments of inertia of the triaxial satellite are denoted by A, B and C, satisfying A ≤ B ≤ C. For convenience, we take the orbital semimajor axis a and the mass of the central planet as the units of length and mass, respectively, and take the time unit to make the mean motion n be unitary.In normalized units, the orbital period of satellite is equal to 2π and the gravitational constant G becomes unitary.
The normalized Hamiltonian, governing the evolution of satellite's spin axis, can be written as (Goldreich & Peale 1966;Flynn & Saha 2005;Gkolias et al. 2016Gkolias et al. , 2019;;Celletti 1990a,b) where θ is the angle between the satellite's longest axis and a reference line, f is the true anomaly and r is the radial distance from the planet and the asphericity parameter α is characterized in terms of the moments of inertia (or equatorial elongation of satellite) as follows (Murray & Dermott 1999): where a s and b s are the semi-axis of triaxial satellite (α is independent upon the shortest semi-major axis c s ).As r(t) and f (t) are functions of time t, the Hamiltonian (1) determines a non-autonomous dynamical model.To make it be autonomous, it is standard to augment the phase space by introducing an action variable T which is conjugated to time t.The resulting new set of canonical variables are denoted by (Flynn & Saha 2005) and the augmented Hamiltonian becomes which determines a 2-DOF dynamical model.In practice, it is usual to take a certain value of p 2 to make the Hamiltonian H(q 1 , p 1 , q 2 , p 2 ) be equal to zero.According to elliptic expansions, the Hamiltonian (2) can be expanded in Fourier series (Celletti 1990a,b).As a consequence, the spin-orbit Hamiltonian truncated at the fourth order in e takes the form: where the coefficients are related to eccentricity in the following form: The equations of motion can be derived from Hamiltonian canonical relations (Morbidelli 2002), In summary, the spin-orbit problem is ruled by two parameters: the orbital eccentricity e and the asphericity parameter α.In Sect.3, the equations of motion are numerically integrated in order to produce dynamical maps as well as Poincaré sections.

DYNAMICAL MAPS
In this section, we numerically explore phase-space structures of the spin-orbit problem by means of dynamical maps, which provides a global view about the dynamics.In particular, two types of indicators are adopted.In producing two types of dynamical maps, the asphericity parameter α and the spin frequency θ are distributed in the interval [0.2, 1.5]×[−0.5,2.5] with steps of δα = 0.01 and δ θ = 0.01 and numerical integration is performed over 500 orbital periods.The first one is the fast Lyapunov indicator (FLI), introduced by Froeschlé et al. (1997), which provides a quick method for distinguishing between chaotic and regular regions in the phase space.There are several variants for the definition of FLI (Barrio 2005(Barrio , 2006;;Skokos 2009).In this work, we adopt the definition of FLI, given by (Guzzo et al. 2002) where X = (θ, θ) T is the state vector, δX (t 0 ) is the initial tangent vector, δX (t) is the deviation at the moment t, and t f is the maximum integration time.In general, the magnitude of FLI increases linearly with time for regular motions and it grows exponentially with time for chaotic motions (Guzzo et al. 2002;Skokos 2009).As a result, the maps of FLI can clearly show regular and chaotic regions.
For the spin-orbit problem, we take the maximum integration time as 100 times the orbital period and fix the initial angle at θ = 0.In practical computation, we take the initial deviation relative to the reference orbit as δX (t 0 ) = √ 2/2 × 10 −4 which include a small component along the tangent vector (our experient shows that the choice of δX (t 0 ) has no influence upon the maps of FLI).Fig. 1 shows the maps of FLI in the ( θ, α) space for two configurations with orbital eccentricities at e = 0.01 and e = 0.02.The blue color, standing for low magnitude of FLI, indicates regular dynamics and the yellow color, corresponding to high magnitude of FLI, indicates chaotic dynamics.From Fig. 1, we can see that (a) the structures corresponding to e = 0.01 and e = 0.02 are qualitatively similar and the latter one has wider chaotic layers than the former one, (b) the chaotic strips in both panels presents V-shape structure in the ( θ, α) space, (c) inside the V-shape region some Lei sub-structures caused by secondary resonances can be visibly observed and (d) outside the V-shape region additional sub-structures caused by high-order spin-orbit resonances are visible.As shown in the following sections, the V-shape strips arising in the maps of FLI are due to the dynamical separatrix of the synchronous primary resonance.Around the separatrix, chaotic layers form the main V-shape structures.
The second indicator we adopt is the secondderivative-based index ∥∆D∥ where D measures the diameter of an orbit, developed by Daquin & Charalambous (2023).The second-derivative-based index is a sensitive and robust indicator, which can be used to detect separatrices, resonant webs and chaotic seas in the phase space.In particular, such an indicator has a strong ability to distinguish minute structures in the phase space.
For the spin-orbit problem, the action-angles variables are denoted by (θ, θ) and the diameter metric of an orbit is defined by where (θ 0 , θ0 ) stand for the initial state and t 0 and t f are the initial and final moments of time.Based on the distance metric D, a normalized second-order-derivativebased quantity is introduced by It should be mentioned that the introduction of the index ∥∆D∥ follows from recent developments on Lagrangian descriptors and arc-length of orbits (Daquin et al. 2022).Regarding the computation method (central difference) about the second derivatives of the diameter, please refer to Daquin & Charalambous (2023) for more details.Contrarily to the FLI method, the method based on ∥∆D∥ is free of tangent dynamics.
The ∥∆D∥ metric allows to sharply detect the separatrix crossing and thus the maps of ∥∆D∥ are able to provide more details on minute structures in the phase space.The difference of ∥∆D∥ defined in this work from the conventional definition given in Daquin & Charalambous (2023) lies in the normalization.Our experiments show that the normalized indicator can detect those structures with small strength.
For the spin-orbit problem, we take the maximum integration time as 100 times the orbital period and fix the initial angle at θ 0 = 0 (this setting is the same as that made in FLI maps).Fig. 2 presents the dynamical structures of the normalized second-derivativebased index ||∆D|| in the ( θ, α) space for eccentricities at e = 0.01 and e = 0.02.Besides the V-shape structures, some minute structures can be visibly detected.
In particular, the structures on the left side of the Vshape structure can be visibly observed in the maps of ||∆D||, while they are not visible in the maps of FLI.In addition, more sub-structures can be found inside the V-shape region.
To understand the structures arising in the dynamical maps shown in Figs. 1 and 2, we take two examples with α = 0.35 and α = 0.65, which are explicitly marked in the left panels.To produce Poincaré sections, we record the states every time when the variable q 2 (or t) satisfies mod (q 2 , 2π) = mod (t, 2π) = 0.
In simulation, the state variable p 2 is chosen to ensure that H is equal to zero and the numerical integration is taken over 500 orbital periods.Poincaré sections of these two representative systems are produced in the case of e = 0.01, as shown in Fig. 3. From the Poincaré sections, we can see that (a) in the case of α = 0.35, the secondary 3:1 and 4:1 spin-orbit resonances can be observed and (b) in the case of α = 0.65 only the secondary 2:1 resonance can be found inside the synchronous primary resonance.These secondary resonance are responsible for the structures arising inside the V-shape region shown in Figs. 1 and 2. In both panels of Fig. 3, we can further observe (a) chaotic layers around the separatrices of the primary resonance, corresponding to the V-shape main structures, and (b) high-order resonances outside the primary resonance, corresponding to the sub-structures outside the V-shape region.
From the dynamical maps shown in Figs. 1 and 2, we may ask: how do high-order and/or secondary spinorbit resonances sculpt the dynamical structures arising in the phase space?In the coming section, we aim at addressing this problem by taking advantage of the theory of perturbative treatments developed by Henrard & Lemaitre (1986) and Henrard (1990).

PERTURBATIVE TREATMENTS
To describe the synchronous primary resonance, we introduce the following canonical transformation: where the argument σ 1 is the resonant angle of the synchronous primary resonance.Under such a new set of variables, the Hamiltonian can be re-organized as fol- lows: From the viewpoint of perturbative treatment (Henrard & Lemaitre 1986;Henrard 1990), the Hamiltonian is composed of the kernel Hamiltonian (unperturbed part) and the part of perturbation, denoted by where the unperturbed part is and the perturbation part is The reason why we keep those periodic terms related to σ 1 in the kernel Hamiltonian is that the magnitudes of  14).The eccentricity is taken as e = 0.01.The red dashed lines stand for the dynamical separatrix of the synchronous primary resonance, specified by equation ( 10).The primary resonance happens inside the region bounded by the dynamical separatrices.
coefficients related to σ 1 are at the zeroth order in e.In this way, the influence of synchronous primary resonance upon high-order and/or secondary resonances are taken into consideration.
The dynamical model specified by the kernel Hamiltonian H 0 is of one degree of freedom and it can be used to describe the dynamical structure of the synchronous primary resonance.Under the perturbation of H 1 , high-order spin-orbit resonances (outside the primary resonance) as well as secondary resonances (inside the primary resonance) may appear in the phase space.
In this section, we will first study the basic dynamics governed by the kernel Hamiltonian and determine the nominal location of high-order and secondary resonances in Sect.4.1, then formulate resonant Hamiltonian models by taking advantage of perturbative treatments (Henrard & Lemaitre 1986;Henrard 1990) in Sect.4.2 and results are shown in Sect.4.3.

Dynamics under the Kernel Hamiltonian model
As σ 2 is absent from H 0 , the unperturbed Hamiltonian model holds an integral of motion Σ 2 .Thus, the kernel Hamiltonian determines an integrable dynamical model.It is not difficult to get that the stable equilibrium points are located at (2σ 1 = 0, Σ 1 = 1) and the unstable equilibrium points are located at (2σ 1 = π, Σ 1 = 1).Thus, the synchronous primary resonance is centered at (2σ 1 = 0, Σ 1 = 1).The level curve of Hamiltonian passing through the saddle point in the phase space, satisfying As a consequence, the dynamical separatrix of the synchronous primary resonance can be explicitly expressed by where the symbol '+' corresponds to the upper separatrix and '-' corresponds to the lower one.Usually, the distance between the upper and lower separatrices evaluated at the resonance centre (2σ 1 = 0, Σ 1 = 1) is called the resonant width.According to equation ( 10), the resonant width (in normalised units) is equal to It shows that the asphericity parameter α stands for the resonant width (or strength) of the synchronous primary resonance.Dynamical separatrix plays an role in dividing the entire phase space into regions of circulation and libration.In particular, the region bounded by the separatrices described by equation ( 10) are referred to as the librating zone and the remaining regions are called circulating zones.In addition, the level curves of Hamiltonian in the phase space (i.e., phase portraits) stand for the trajectories, which are exactly periodic.
Under the unperturbed Hamiltonian model, let us introduce action-angle variables as follows (Morbidelli 2002): where σ * 1 becomes a linear function of time, ) is a periodic function with period T and Σ * 1 is called Arnold action, which represents the signed area bounded by the level curves in the phase space divided by 2π (Morbidelli 2002).The Arnold action is equal to zero at the centre of the synchronous resonance and it reaches the maximum along the separatrix of the primary resonance.
The transformation given by equation ( 11) is canonical with the following generating function, After the canonical transformation, the unperturbed Hamiltonian can be written as a normal form, which determines the fundamental frequencies, given by In particular, when the fundamental frequencies become commensurable, high-order or secondary spin-orbit resonances may happen.As a result, the condition of a certain k 1 :k 2 spin-orbit resonance takes the form, By solving equation ( 14), it is possible to determine the nominal location of high-order and/or secondary resonances.
In Fig. 3, the nominal locations of high-order and secondary spin-orbit resonances are shown in the phase space (θ, θ) for the spin-orbit problems with α = 0.35 and α = 0.65 and they are compared with the associated Poincaré sections.In both panels of Fig. 3, the red dashed lines stand for the dynamical separatrices of the synchronous primary resonance, which are mathematically described by equation (10).It is observed that around the dynamical separatrix there is a chaotic layer, which is wilder with a higher asphericity parameter α.
For convenience, we refer to the spin-orbit resonances inside the primary resonance as the secondary resonances and the ones outside the primary resonance as high-order resonances.In the case of α = 0.35, the secondary 3:1 and 4:1 resonances can be observed and, in the case of α = 0.65, only the secondary 2:1 resonance can be seen.Comparing the distribution of high-order and secondary resonances to the Poincaré sections, we can further observe that the distributions of high-order and secondary resonances are in good agreement with the structures arising in Poincaré sections.
In Fig. 4, we further present the distribution of highorder and secondary resonance in the ( θ, α) space and compare them to the dynamical map of the normalized second-derivative-based index ||∆D||.The map of ||∆D|| shown in the right panel is the same as the one shown in Fig. 2. In the left panel, the red shaded region stands for the space where the synchronous primary resonance takes place and the red dashed lines stand for the separatrices of the primary resonance.Secondary resonances happen inside the shaded region and highorder resonances takes place outside the shaded region.It is observed that the nominal locations of spin-orbit resonances are symmetric with respect to the nominal location of the primary resonance at θ = 1.
From the right panel of Fig. 4, we can observe that the distributions of high-order and secondary resonances determined under the unperturbed Hamiltonian model are in excellent agreement with the structures arising in the map of ||∆D||.

Resonant Hamiltonian model
From the previous subsection, it is known that highorder and secondary resonances may take place in the phase space.The purpose of this section is to formulate a Hamiltonian model for describing the resonant dynamics.
Under the transformation given by equation ( 11), the Hamiltonian can be re-organized as where the perturbation part H 1 is much smaller than the Kernel Hamiltonian H 0 .Therefore, it is reasonable to take H 1 as a perturbation.
It should be mentioned that the core concept of henrard's perturbative treatment (Henrard & Lemaitre 1986;Henrard 1990) lies in transforming the Hamiltonian form given by equation ( 9) to the standard form of Hamiltonian given by equation ( 15).In the standard Hamiltonian, the kernel Hamiltonian is of normal form.The standard Hamiltonian form given by equation ( 15) is ready for studying resonances between σ * 1 and σ * 2 .In order to study the k 1 :k 2 resonance between σ * 1 and σ * 2 , let us introduce the following transformation, which is canonical with the generating function The resonant argument for the Under the new set of canonical variables given by equation ( 16), the Hamiltonian can be further expressed as

Lei
In particular, when the triaxial satellite is located inside a certain k 1 :k 2 spin-orbit resonance, the corresponding resonant argument γ 1 = σ * 1 − k2 k1 σ * 2 becomes a longperiodic angle and the other one γ 2 is a short-periodic angle.It means that equation ( 17) determines a separable Hamiltonian system (Henrard 1990).As a consequence, in the long-term evolution, those short-periodic effects can be filtered out by means of averaging theory (corresponding to the lowest-order perturbation theory).To this end, we perform average for the Hamiltonian over k 1 times the period of γ 2 , leading to the resonant Hamiltonian, Under the resonant Hamiltonian model determined by H * , the angle γ 2 becomes a cyclic variable, so that its conjugate momentum is a motion integral.As a result, the resonant Hamiltonian model determined by equation ( 18) is integrable.
Given the motion integral, the global dynamical structures can be explored by plotting level curves of resonant Hamiltonian H * in the phase space (i.e., phase portraits).In order to compare phase portraits with Poincaré sections, we need to project the phase space (σ * 1 , Σ * 1 ) to the space (θ, θ) by considering the condition of Poincaré sections: mod (q 2 , 2π) = 0.

Results
Analytical and numerical results for high-order and secondary spin-orbit resonances in the phase space (θ, θ) are presented in Fig. 5. Regarding analytical result, we provide the dynamical separatrix, which stands for the isolines of resonant Hamiltonian.Dynamical separatrix divides the phase space into regions of circulation and libration.In essence, dynamical separatrix corresponds to invariant manifolds stemming from unstable periodic orbits in the original dynamical model.In the left-column panels, dynamical separatrices of the synchronous primary resonance, high-order resonances and secondary resonances are shown and, in the right-column panels, they structures are compared to numerical structures arising in Poincaré sections.In practical simulations, the spin-orbit problems with asphericity parameters at α = 0.35, α = 0.65, α = 1.1 and α = 1.4 are considered and the eccentricity is taken as e = 0.01.In all panels, the dynamical separatrices of the synchronous primary resonance are marked in red dashed lines.
From Fig. 5, it is observed that the analytical structures (i.e., dynamical separatrices) of high-order and/or secondary spin-orbit resonances are in good agreement with numerical structures arising in the Poincaré sections.In particular, the numerical structures are dominated by the secondary 3:1 and 4:1 spin-orbit resonances and high-order resonances including 2:1, 4:1 and 2:−1 in the case of α = 0.35, by the secondary 2:1 resonance and high-order resonances including 2:1 and 2:−1 in the case of α = 0.65 and by the secondary 1:1 resonance and high-order 1:1 resonance in the cases of α = 1.1 and α = 1.4.
In the spin-orbit problem with α = 0.35, the dynamical separatrices of the secondary 4:1 resonance (inside the primary resonance) and high-order 4:1 resonance (outside the primary resonance) are close to the separatrices of the primary resonance, and they provide boundaries for the chaotic layer around the separatrix of the primary resonance.When the asphericity parameter is increased to α = 0.65, the high-order 4:1 resonance disappears and the secondary 3:1 and 4:1 resonances are replaced by the secondary 2:1 resonance.It is observed that the dynamical separatrices of the high-order 2:1 and secondary 2:1 resonances provide boundaries for chaotic layer around the separatrix of the primary resonance.When the asphericity parameter is further increased up to α = 1.1 and α = 1.4,only high-order 1:1 and secondary 1:1 resonances can be observed in the phase space.
According to the resonant Hamiltonian model, we can analytically determine the location of synchronous primary and/or secondary 1:1 resonances (characteristic curves).On the other hand, the numerical characteristic curves can be determined by computing the associated periodic orbits in a high accuracy.Denote the period by T .Periodic orbits should satisfy the following condition: Without loss of generality, we assume the angle at t = t 0 is equal to zero, i.e. θ 0 = 0. Thus, there are only two variables ( θ0 , T ) to be determined for a periodic orbit.In practice, we adopt the Newton-Raphson method to solve the unknown variables of periodic orbits.
Fig. 6 shows the analytical characteristic curves in the left panel and comparisons between analytical and numerical curves in the right panel.The shaded region stands for the space where the synchronous primary resonance happens.From Fig. 6, we can see that (a) the analytical characteristic curves can agree well with the numerical ones in the entire considered range α ∈ [0, 1.5], and (b) the secondary 1:1 resonance appears when the asphericity parameter is higher than α∼1.06 in the case of e = 0.01.
Regarding spin-orbit problem, Gkolias et al. (2016Gkolias et al. ( , 2019) ) took advantage of Lie-series transformation theory to study the low-order secondary resonances bifurcating from the synchronous primary resonance.Analytical characteristic curves are produced from the 11th normal form construction as well as from wisdom formula derived from nonlinear method of Bogoliubov and Mitropolsky (Wisdom et al. 1984), and analytical results are compared to numerical computations.It is concluded that the normal form method has a significant improvement compared to the nonlinear method of Bogoliubov and Mitropolsky.In particular, the normal form solution can agree with numerical result in the range of α ∈ [0, 1.2], but starts to diverge when α is greater than 1.2 (please refer to Fig. 2 in their work for more details).However, our analytical estimate can match the numerical solution in the entire considered range of α up to 1.5 (see the right panel of Fig. 6).
Fig. 7 shows all the analytical characteristic curves of the synchronous primary resonance, secondary resonances and high-order resonances in the left panel and makes a comparison between analytical curves to the numerical map of ||∆D|| in the right panel.We can see that the structures arising in the dynamical map have a good correspondence to the analytical characteristic curves.Thus, these numerical structures can be well understood thanks to the associated spin-orbit resonance dynamics.
In particular, Fig. 7 provides a global view for the rotational dynamics of the secondary in the ( θ, α) space.It is possible for us to predict possible dynamical outcomes for binary asteroid systems according to the associated asphericity parameters.For example, the following predictions can be made: (a) secondary 1:1 resonance may appear when α > 1.06, (b) high-order 2:1 and 2:−1 resonances as well as secondary 2:1, 3:1, and 4:1 resonances may disappear in the high-α space because they are inside the chaotic sea around the separatrix of the synchronous primary resonance, and (c) high-order 1:1 and 1:−1 resonances can appear in the entire considered space but the strength of high-order 1:−1 resonance is relatively weak.

APPLICATIONS TO BINARY ASTEROID SYSTEMS
The analytical approach is applied to two representative binary asteroid systems: (65803) Didymos and ( 4383) Suruga (Pravec et al. 2016;Warner 2013).( 65803) Didymos is the target of the DART mission (Agrusa et al. 2022).It should be mentioned that the analytical method presented in this work is also applicable for other binary asteroid systems.The loca- .Analytical characteristic curves of the synchronous primary resonance and the secondary 1:1 spin-orbit resonance (left panel ) and comparisons to the associated numerical curves (right panel ).The vertical dashed line stands for the nominal location of the synchronous primary resonance and the shaded area represents the space where the synchronous primary resonance takes place.The agreement between analytical and numerical curves is good over the entire considered range α ∈ [0, 1.5].
Table 1.Physical parameters for the example binary asteroid systems adopted in this work (please refer to http://www.asu.cas.cz/asteroid/binastdata.htm for detailed data of binary asteroid systems).The ratio of spin angular momentum to the orbital angular momentum for the secondary Lspin/L orb is estimated from the method given in Jafari- Nadoushan (2023)  tion of these two systems are marked in the right panel of Fig. 7, which shows that their asphericity parameters are greater than the critical value α c = 1.06 at which the secondary 1:1 resonance bifurcates from the synchronous primary resonance.The physical parameters of the adopted binary asteroid systems are shown in Table 1, including the ratio of the secondary diameter to the primary diameter D s /D p , ratio of the secondary mass to the primary mass m s /m p , spin period of the primary P p , mutual orbital period of system P orb , spin period of the secondary P s , the equatorial elongation of the secondary a s /b s , ratio of the orbital semimajor axis to the secondary's diameter a orb /D s , mutual orbital eccentricity e orb , ratio of spin angular momentum to the orbital angular momentum L spin /L orb and asphericity parameter α.
Regarding binary asteroid systems, Jafari-Nadoushan (2023) offered a criterion for the definition of dynamical closeness of asteroid pairs: a binary asteroid system is a close pair if the ratio of rotational angular momentum to the orbital angular momentum is greater than 10 percent.According to their estimate, we can get the ratios of spin angular momentum to the orbital angular momentum are L spin /L orb = 2.2 × 10 −3 for (65803) Didymos and L spin /L orb = 9.9 × 10 −4 for (4383) Suruga, which are much smaller than 10 percent.Thus, both binary asteroid systems are far away from close binaries, which means that the couple between orbital and rotational evolution is very weak.In the dynamical evolution, the influence of rotational evolution upon the orbital dynamics can be ignored.As a result, the spin-orbit problem discussed in this work is a good approximation for such a kind of weakly-coupled binary systems.
According to Table 1, the asphericity parameter is α = 1.1 for binary asteroid system (65803) Didymos and it is equal to α = 1.36 for binary asteroid system (4383) Suruga.Considering that the orbital eccentricities are small and highly uncertain, we take two small eccentricities at e = 0.01 and e = 0.02 in practical simulations.
Analytical and numerical results are presented in Fig. 8, where the top panels are for (65803) Didymos and the bottom panels are for (4383) Suruga.The dynamical separatrix of the synchronous primary resonance is shown in red dashed line and the separatrix of the secondary 1:1 resonance is given in blue line.From Fig. 8, we can observe that (a) there is an excellent agreement between the analytical and numerical structures, and (b) chaotic layer becomes wilder with a higher eccentricity.
According to the current physical parameters, it is highly possible that the binary asteroid systems (65803) Didymos and (4383) Suruga are inside the secondary 1:1 spin-orbit resonance.However, there are a large uncertainty for their elongation a s /b s , leading to a large uncertainty for determining asphericity parameter α.For example, it will be outside the secondary 1:1 resonance if the uncertainty of asphericity parameter makes α be smaller than ∼1.06.Thus, more observations are required to further determine the asphericity parameter α accurately in order to predict their real resonant states.However, Fig. 7 can help to predict other possible resonant states when α is varied.

CONCLUSIONS
In this work, high-order and secondary resonances in the spin-orbit problem are investigated.Dynamical maps are numerically produced in order to provide a global picture for the dynamics.In particular, two types of dynamical indicators are introduced.The first one is the fast Lyapunov indicator (FLI) and the second one is the normalized second-derivative-based index ||∆D||.The maps of FLI can distinguish the regular and chaotic regions and the maps of ||∆D|| has an ability to detect minute structures in the phase space.It is concluded that the main V-shape structures are sculpted by the synchronous primary resonance, those minute structures inside the V-shape region are governed by secondary resonances and the ones outside the V-shape region are dominated by high-order spin-orbit resonances.
The main purpose of this work is to understand dynamical causes of numerical structures.To this end, the theory of perturbative treatments is adopted.In particular, the Hamiltonian is divided into the kernel Hamiltonian and the perturbation part.The terms related to the synchronous primary resonance are included in the kernel Hamiltonian, so that it is possible to consider the influence of the synchronous primary resonance upon high-order and/or secondary resonances.Under the kernel Hamiltonian model, a series of canonical transformations are introduced and nominal locations of high-order and secondary resonances are determined.By taking advantage of Henrard's perturbative treatments (Henrard 1990), resonant Hamiltonian model is formulated.Thus, it is ready to study high-order and secondary resonance dynamics analytically.Results show that there is an excellent correspondence between analytical structures arising in phase portraits and numerical structures arising both in the dynamical maps and in the Poincaré sections.Applications to two binary asteroid systems (including Didymos and Suruga) show that these two binary systems are possible to be inside the secondary 1:1 spin-orbit resonance.4383) Suruga (bottom panels).The blue lines stand for the dynamical separatrices of the 1:1 secondary resonance and red dashed lines stand for the dynamical separatrices of the synchronous primary resonance.The asphericity parameter is α = 1.1 for binary asteroid system (65803) Didymos and it is α = 1.36 for binary asteroid system (4383) Suruga.
The main results are presented in Fig. 7 where the distributions of synchronous primary resonance, high-order resonance and secondary resonances are presented and compared with dynamical map of ||∆D||.It provides a global picture for us to predict possible resonant states when the asphericity parameter is changed.

Figure 1 .Figure 2 .
Figure 1.Dynamical maps of the fast Lyapunov indicator (FLI) shown in the ( θ, α) space for the configurations with orbital eccentricities at e = 0.01 (left panel ) and e = 0.02 (right panel ).The index shown in the color bar corresponds to the base 10 logarithm of FLI.In both panels, the main V-shape structures are visible.The red dashed lines stand for the locations of α = 0.35 and α = 0.65 and their associated Poincaré sections are shown in Fig. 3.

Figure 3 .
Figure 3. Nominal location of the centers of high-order and secondary resonances determined under the Kernel Hamiltonian model together with the Poincaré surfaces of section for spin-orbit problems with asphericity parameters at α = 0.35 (left panel ) and α = 0.65 (right panel ).Nominal locations of high-order and secondary resonances are determined by equation (14).The eccentricity is taken as e = 0.01.The red dashed lines stand for the dynamical separatrix of the synchronous primary resonance, specified by equation (10).The primary resonance happens inside the region bounded by the dynamical separatrices.

Figure 4 .
Figure 4. Nominal locations of high-order and secondary spin-orbit resonances determined under the Kernel Hamiltonian model shown in the ( θ, α) space (left panel ) and comparisons to dynamical maps of the second-derivative-based index in the case of e = 0.01 (right panel ).The map shown in the right panel is the same as that presented in the left panel of Fig. 2. In the left panel, the shaded region stands for the parameter space where the synchronous primary resonance happens.Those resonances inside the primary resonance are called secondary resonances and the ones outside the primary resonance are referred to as high-order spin-orbit resonances.

Figure 5 .
Figure5.Dynamical separatrices of high-order and secondary spin-orbit resonances in the phase space (left-column panels) as well as their comparisons to numerical structures arising in Poincaré sections (right-column panels) for spin-orbit problems with different asphericity parameters.The red dashed lines stand for dynamical separatrices of the synchronous primary resonance, which are expressed by equation (10).
Figure6.Analytical characteristic curves of the synchronous primary resonance and the secondary 1:1 spin-orbit resonance (left panel ) and comparisons to the associated numerical curves (right panel ).The vertical dashed line stands for the nominal location of the synchronous primary resonance and the shaded area represents the space where the synchronous primary resonance takes place.The agreement between analytical and numerical curves is good over the entire considered range α ∈ [0, 1.5].

Figure 7 .
Figure7.Analytical characteristic curves of the synchronous primary resonance, high-order resonances (outside the primary resonance) and secondary resonances (inside the primary resonance) in the ( θ, α) space (left panel ) and comparisons with dynamical maps of the normalized second-derivative-based index ||∆D|| (right panel ).In the right panel, the locations of two typical binary-asteroid systems are marked.The asphericity parameter is α = 1.1 for the binary-asteroid system (65803) Didymos and it is 1.36 for (4383) Suruga.The vertical dashed line stands for the nominal location of the synchronous primary resonance and the shaded region in the left panel represents the space where the synchronous primary resonance happens.

Figure 8 .
Figure8.Analytical structures (dynamical separatrices) of the synchronous primary resonance and the secondary 1:1 spinorbit resonance together with the associated numerical structures arising in Poincaré sections at different eccentricities for binary asteroid systems (65803) Didymos (top panels) and (4383) Suruga (bottom panels).The blue lines stand for the dynamical separatrices of the 1:1 secondary resonance and red dashed lines stand for the dynamical separatrices of the synchronous primary resonance.The asphericity parameter is α = 1.1 for binary asteroid system (65803) Didymos and it is α = 1.36 for binary asteroid system (4383) Suruga.
This work is supported by the National Natural Science Foundation of China (Nos.12073011, and 12233003) and the National Key R&D Program of China (No. 2019YFA0706601). .