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

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 (Peale & Gold 1965; Lemaitre et al. 2006) 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 (Scheeres et al. 2006; Naidu & Margot 2015; Pravec et al. 2016, 2019).

The conventional model of spin–orbit resonance is formulated under the assumption that the two bodies move around their barycenter in fixed Keplerian orbits (Goldreich & Peale 1966; Celletti 1990a, 1990b). 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

$$\begin{eqnarray*}&&\alpha \gt {\alpha }_{c}=\displaystyle \frac{1}{2+\sqrt{14e}}.\end{eqnarray*}$$

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 structures by taking advantage of Poincaré sections. Later, Nadoushan & Assadian (2016) provided a geography of varieties of rotational resonances under a fourth-order Hamiltonian model.

If the orbit and rotation of a binary are coupled, then it becomes a 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 spin–spin 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 an 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 the primary spin–orbit resonances. However, if the frequency of small-amplitude oscillations about the primary resonance and the orbital frequency are commensurable, then secondary spin–orbit resonances may happen. The existence of secondary resonances make the phase-space structures 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. (2016, 2019) took advantage of Lie-series transformation theory to formulate analytical approximations for those low-order 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 the 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. (2016, 2019) can be used to study secondary resonances but cannot be used to study those high-order spin–orbit resonances. Thus, a unified resonant Hamiltonian model for describing high-order and secondary resonances in the spin–orbit problem is required to explain phase-space structures.

In this work, two types of dynamical maps (i.e., the fast Lyapunov indicator and the normalized second-derivative-based index) are produced in the $(\dot{\theta },\alpha )$ 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.

The remaining part of this work is organized as follows. In Section 2, the Hamiltonian function of the spin–orbit 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 the conclusions are summarized in Section 6.

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 spin axis (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 ${ \mathcal G }$ becomes unitary.

The normalized Hamiltonian, governing the evolution of satellite's spin axis, can be written as (Goldreich & Peale 1966; Celletti 1990a, 1990b; Flynn & Saha 2005; Gkolias et al. 2016, 2019)

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \frac{1}{2}{\dot{\theta }}^{2}-\displaystyle \frac{{\alpha }^{2}}{4{r}^{3}}\cos \left[2\left(\theta -f\right)\right],\end{eqnarray} \tag{ 1 }$$

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. The asphericity parameter α is characterized in terms of the moments of inertia (or equatorial elongation of satellite) as follows (Murray & Dermott 1999):

$$\begin{eqnarray*}&&\alpha =\sqrt{\displaystyle \frac{3\left(B-A\right)}{C}}=\sqrt{\displaystyle \frac{3\left({a}_{s}^{2}/{b}_{s}^{2}-1\right)}{{a}_{s}^{2}/{b}_{s}^{2}+1}},\end{eqnarray*}$$

where a_s and b_s are the semi-axis of triaxial satellite (α is independent upon the shortest semimajor axis c_s ). Since 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)

$$\begin{eqnarray*}&&\begin{array}{l}{q}_{1}=\theta ,\quad {p}_{1}=\dot{\theta },\\ {q}_{2}=t,\quad {p}_{2}=T\end{array}\end{eqnarray*}$$

and the augmented Hamiltonian becomes

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \frac{1}{2}{p}_{1}^{2}+{p}_{2}-\displaystyle \frac{{\alpha }^{2}}{4{r}^{3}}\cos \left[2\left({q}_{1}-f\right)\right],\end{eqnarray} \tag{ 2 }$$

which determines a 2-DOF dynamical model. In practice, it is usual to take a certain value of p₂ to make the Hamiltonian ${ \mathcal H }({q}_{1},{p}_{1},{q}_{2},{p}_{2})$ be equal to zero.

According to elliptic expansions, the Hamiltonian (2) can be expanded as a Fourier series (Celletti 1990a, 1990b). As a consequence, the spin–orbit Hamiltonian truncated at the fourth order in e takes the form:

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \frac{{p}_{1}^{2}}{2}+{p}_{2}-\displaystyle \frac{{\alpha }^{2}}{4}\displaystyle \sum _{n=-6}^{2}{{ \mathcal C }}_{n}(e)\cos \left(2{q}_{1}+{{nq}}_{2}\right),\end{eqnarray} \tag{ 3 }$$

where the coefficients are related to eccentricity in the following form:

$$\begin{eqnarray*}&&\begin{array}{l}{{ \mathcal C }}_{-6}=\displaystyle \frac{533}{16}{e}^{4},\ {{ \mathcal C }}_{-5}=\displaystyle \frac{845}{48}{e}^{3},\ {{ \mathcal C }}_{-4}=\displaystyle \frac{17}{2}{e}^{2}-\displaystyle \frac{115}{6}{e}^{4},\\ {{ \mathcal C }}_{-3}=\displaystyle \frac{7}{2}e-\displaystyle \frac{123}{16}{e}^{3},\ {{ \mathcal C }}_{-2}=1-\displaystyle \frac{5}{2}{e}^{2}+\displaystyle \frac{13}{16}{e}^{4},\\ {{ \mathcal C }}_{-1}=-\displaystyle \frac{1}{2}e+\displaystyle \frac{1}{16}{e}^{3},\ {{ \mathcal C }}_{0}=0,\ {{ \mathcal C }}_{1}=\displaystyle \frac{1}{48}{e}^{3},\ {{ \mathcal C }}_{2}=\displaystyle \frac{1}{24}{e}^{4}.\end{array}\end{eqnarray*}$$

The equations of motion can be derived from Hamiltonian canonical relations (Morbidelli 2002),

$$\begin{eqnarray}&&\begin{array}{l}{\dot{q}}_{1}=\displaystyle \frac{\partial { \mathcal H }}{\partial {p}_{1}},\quad {\dot{p}}_{1}=-\displaystyle \frac{\partial { \mathcal H }}{\partial {q}_{1}},\\ {\dot{q}}_{2}=\displaystyle \frac{\partial { \mathcal H }}{\partial {p}_{2}},\quad {\dot{p}}_{2}=-\displaystyle \frac{\partial { \mathcal H }}{\partial {q}_{2}}.\end{array}\end{eqnarray} \tag{ 4 }$$

In summary, the spin–orbit problem is ruled by two parameters: the orbital eccentricity e and the asphericity parameter α. In Section 3, the equations of motion are numerically integrated to produce dynamical maps as well as Poincaré sections.

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 $\dot{\theta }$ are distributed in the interval [0.2, 1.5] × [−0.5, 2.5] with steps of δ α = 0.01 and $\delta \dot{\theta }=0.01$, and numerical integration is performed over 500 orbital periods.

The first 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, 2006; Skokos 2009). In this work, we adopt the definition of FLI, given by Guzzo et al. (2002)

$$\begin{eqnarray}&&\mathrm{FLI}=\mathop{\sup }\limits_{t\in \left[{t}_{0},{t}_{f}\right]}{\mathrm{log}}_{10}\displaystyle \frac{\parallel \delta {\boldsymbol{X}}\left(t\right)\parallel }{\parallel \delta {\boldsymbol{X}}\left({t}_{0}\right)\parallel }\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{X}}={(\theta ,\dot{\theta })}^{{\rm{T}}}$ is the state vector, $\delta {\boldsymbol{X}}\left({t}_{0}\right)$ is the initial tangent vector, $\delta {\boldsymbol{X}}\left(t\right)$ is the deviation at the moment t, and t_f is the maximum integration time. In general, the norm of a deviation vector increases linearly with time for regular motions and it grows exponentially with time for chaotic motions. As a result, it is possible to discriminate regular and chaotic orbits according to the magnitude of FLI (Guzzo et al. 2002; Skokos 2009).

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 norm of initial deviation relative to the reference orbit as $| | \delta {\boldsymbol{X}}\left({t}_{0}\right)| | \,=\sqrt{2}/2\,\times \,{10}^{-4}$, which includes a small component along the tangent vector (our experiment shows that the choice of $\delta {\boldsymbol{X}}\left({t}_{0}\right)$ has no influence upon the maps of FLI). Figure 1 shows the maps of FLI in the $(\dot{\theta },\alpha )$ 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 Figure 1, we can see that (a) the structures corresponding to e = 0.01 and e = 0.02 are qualitatively similar and the latter has wider chaotic layers than the former, (b) the chaotic strips in both panels presents V-shape structure in the $(\dot{\theta },\alpha )$ space, (c) inside the V-shape region some substructures caused by secondary resonances can be visibly observed, and (d) outside the V-shape region additional substructures 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.

Zoom In Zoom Out Reset image size

The second indicator that we adopt is the second-derivative-based index $\parallel {\rm{\Delta }}D\parallel$ where D measures the diameter of an orbit, which was 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–angle variables are denoted by $(\theta ,\dot{\theta })$ and the diameter metric of an orbit is defined by

$$\begin{eqnarray}&&D\left({\theta }_{0},{\dot{\theta }}_{0}\right)=\mathop{\max }\limits_{{t}_{0}\leqslant t\leqslant {t}_{f}}\dot{\theta }\left({\theta }_{0},{\dot{\theta }}_{0},t\right)-\mathop{\min }\limits_{{t}_{0}\leqslant t\leqslant {t}_{f}}\dot{\theta }\left({\theta }_{0},{\dot{\theta }}_{0},t\right)\end{eqnarray} \tag{ 6 }$$

where $({\theta }_{0},{\dot{\theta }}_{0})$ stands for the initial state, and t₀ and t_f are the initial and final moments of time. Based on the distance metric D, a normalized second-order-derivative-based quantity is introduced by

$$\begin{eqnarray}&&\parallel {\rm{\Delta }}D\parallel =\displaystyle \frac{1}{D\left({\theta }_{0},{\dot{\theta }}_{0}\right)}\left(\left|\displaystyle \frac{{\partial }^{2}D\left({\theta }_{0},{\dot{\theta }}_{0}\right)}{\partial {\dot{\theta }}_{0}^{2}}\right|+\left|\displaystyle \frac{{\partial }^{2}D\left({\theta }_{0},{\dot{\theta }}_{0}\right)}{\partial {\theta }_{0}^{2}}\right|\right).\end{eqnarray} \tag{ 7 }$$

It should be mentioned that the introduction of the index $\parallel {\rm{\Delta }}D\parallel$ 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 $\parallel {\rm{\Delta }}D\parallel$ is free of tangent dynamics.

The $\parallel {\rm{\Delta }}D\parallel$ metric allows to sharply detect the separatrix crossing, and thus the maps of $\parallel {\rm{\Delta }}D\parallel$ are able to provide more details on minute structures in the phase space. The difference of $\parallel {\rm{\Delta }}D\parallel$ 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 (this setting is the same as that made in FLI maps). Figure 2 presents the dynamical structures of the normalized second-derivative-based index ∣∣ΔD∣∣ in the $(\dot{\theta },\alpha )$ 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-hand side of the V-shape structure can be visibly observed in the maps of ∣∣ΔD∣∣, while they are not visible in the maps of FLI. In addition, more substructures can be found inside the V-shape region.

Zoom In Zoom Out Reset image size

To understand the structures arising in the dynamical maps shown in Figures 1 and 2, we take two examples with α = 0.35 and α = 0.65 to produce Poincaré sections by recording the states every time when the variable q₂ (or t) satisfies

$$\begin{eqnarray*}&&\mathrm{mod}\left({q}_{2},2\pi \right)=\mathrm{mod}\left(t,2\pi \right)=0.\end{eqnarray*}$$

In simulation, the state variable p₂ is chosen to ensure that ${ \mathcal 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 Figure 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 resonances are responsible for the structures arising inside the V-shape region shown in Figures 1 and 2. In both panels of Figure 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 substructures outside the V-shape region.

Zoom In Zoom Out Reset image size

From the dynamical maps shown in Figures 1 and 2, we may ask: how do high-order and/or secondary spin–orbit resonances sculpt the dynamical structures arising in the phase space? In the next section, we aim at addressing this problem by taking advantage of the theory of perturbative treatments developed by Henrard & Lemaitre (1986) and Henrard (1990).

To describe the synchronous primary resonance, we introduce the following canonical transformation:

$$\begin{eqnarray*}&&\begin{array}{l}{\sigma }_{1}={q}_{1}-{q}_{2},\quad {{\rm{\Sigma }}}_{1}={p}_{1},\\ {\sigma }_{2}={q}_{2},\quad {{\rm{\Sigma }}}_{2}={p}_{2}+{p}_{1},\end{array}\end{eqnarray*}$$

where the argument σ₁ is the resonant angle of the synchronous primary resonance. Under such a new set of variables, the Hamiltonian can be reorganized as follows:

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal H }=\displaystyle \frac{{{\rm{\Sigma }}}_{1}^{2}}{2}+{{\rm{\Sigma }}}_{2}-{{\rm{\Sigma }}}_{1}-\displaystyle \frac{{\alpha }^{2}}{4}{{ \mathcal C }}_{-2}\left(e\right)\cos \left(2{\sigma }_{1}\right)\\ -\displaystyle \frac{{\alpha }^{2}}{4}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{-6\leqslant n\leqslant 2}{n\ne -2}}{{ \mathcal C }}_{n}\left(e\right)\cos \left[2{\sigma }_{1}+\left(n+2\right){\sigma }_{2}\right].\end{array}\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&{ \mathcal H }={{ \mathcal H }}_{0}\left({\sigma }_{1},{{\rm{\Sigma }}}_{1},{{\rm{\Sigma }}}_{2}\right)+{{ \mathcal H }}_{1}\left({\sigma }_{1},{{\rm{\Sigma }}}_{1},{\sigma }_{2},{{\rm{\Sigma }}}_{2}\right),\end{eqnarray} \tag{ 9 }$$

where the unperturbed part is

$$\begin{eqnarray*}&&{{ \mathcal H }}_{0}=\displaystyle \frac{{{\rm{\Sigma }}}_{1}^{2}}{2}+{{\rm{\Sigma }}}_{2}-{{\rm{\Sigma }}}_{1}-\displaystyle \frac{{\alpha }^{2}}{4}\cos \left(2{\sigma }_{1}\right)\end{eqnarray*}$$

and the perturbation part is

$$\begin{eqnarray*}&&\begin{array}{l}{{ \mathcal H }}_{1}=-\displaystyle \frac{{\alpha }^{2}}{4}\left[{{ \mathcal C }}_{-2}\left(e\right)-1\right]\cos \left(2{\sigma }_{1}\right)\\ -\displaystyle \frac{{\alpha }^{2}}{4}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{-6\leqslant n\leqslant 2}{n\ne -2}}{{ \mathcal C }}_{n}\left(e\right)\cos \left[2{\sigma }_{1}+\left(n+2\right){\sigma }_{2}\right].\end{array}\end{eqnarray*}$$

The reason why we keep these periodic terms related to σ₁ in the kernel Hamiltonian is that the magnitudes of coefficients related to σ₁ 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 ${{ \mathcal 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 ${{ \mathcal 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 Section 4.1. We then formulate resonant Hamiltonian models by taking advantage of perturbative treatments (Henrard & Lemaitre 1986; Henrard 1990) in Section 4.2 and the results are shown in Section 4.3.

4.1. Dynamics Under the Unperturbed Hamiltonian Model

Since σ₂ is absent from ${{ \mathcal H }}_{0}$, the unperturbed Hamiltonian model holds an integral of motion Σ₂. Thus, the kernel Hamiltonian determines an integrable dynamical model. It is not difficult to get that the stable equilibrium points are located at (2σ₁ = 0, Σ₁ = 1) and the unstable equilibrium points are located at (2σ₁ = π, Σ₁ = 1). Thus, the synchronous primary resonance is centered at (2σ₁ = 0, Σ₁ = 1). The level curve of Hamiltonian passing through the saddle point in the phase space, satisfying

$$\begin{eqnarray*}&&\displaystyle \frac{1}{2}\left({{\rm{\Sigma }}}_{1}^{2}+1\right)-{{\rm{\Sigma }}}_{1}-\displaystyle \frac{{\alpha }^{2}}{4}\left(\cos 2{\sigma }_{1}+1\right)=0.\end{eqnarray*}$$

As a consequence, the dynamical separatrix of the synchronous primary resonance can be explicitly expressed by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{1}=1\pm \displaystyle \frac{\sqrt{2}}{2}\alpha \sqrt{1+\cos 2{\sigma }_{1}},\end{eqnarray} \tag{ 10 }$$

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 center (2σ₁ = 0, Σ₁ = 1) is called the resonant width. According to Equation (10), the resonant half width (in normalized units) is equal to

$$\begin{eqnarray*}&&{\rm{\Delta }}{{\rm{\Sigma }}}_{1}={\rm{\Delta }}\dot{\theta }=\alpha .\end{eqnarray*}$$

This shows that the asphericity parameter α stands for the resonant half 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):

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{1}^{* }={\sigma }_{1}-\rho \left({\sigma }_{1}^{* },{{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right),\quad {{\rm{\Sigma }}}_{1}^{* }=\displaystyle \frac{1}{2\pi }\oint {{\rm{\Sigma }}}_{1}{\rm{d}}{\sigma }_{1}\\ {\sigma }_{2}^{* }={\sigma }_{2},\quad {{\rm{\Sigma }}}_{2}^{* }={{\rm{\Sigma }}}_{2},\end{array}\end{eqnarray} \tag{ 11 }$$

where ${\sigma }_{1}^{* }$ becomes a linear function of time, $\rho \left({\sigma }_{1}^{* },{{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right)$ is a periodic function with period T and ${{\rm{\Sigma }}}_{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 center 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,

$$\begin{eqnarray*}&&S\left({\sigma }_{1},{{\rm{\Sigma }}}_{1}^{* },{\sigma }_{2},{{\rm{\Sigma }}}_{2}^{* }\right)={\sigma }_{2}{{\rm{\Sigma }}}_{2}^{* }+\int {{\rm{\Sigma }}}_{1}\left({{ \mathcal H }}_{0}\left({{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right),{\sigma }_{1}\right){\rm{d}}{\sigma }_{1}.\end{eqnarray*}$$

After the canonical transformation, the unperturbed Hamiltonian can be written as a normal form,

$$\begin{eqnarray}&&{{ \mathcal H }}_{0}\left({\sigma }_{1},{{\rm{\Sigma }}}_{1},{{\rm{\Sigma }}}_{2}\right)\Rightarrow {{ \mathcal H }}_{0}\left({{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right),\end{eqnarray} \tag{ 12 }$$

which determines the fundamental frequencies, given by

$$\begin{eqnarray}&&{\dot{\sigma }}_{1}^{* }=\displaystyle \frac{\partial {{ \mathcal H }}_{0}\left({{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right)}{\partial {{\rm{\Sigma }}}_{1}^{* }},\quad {\dot{\sigma }}_{2}^{* }=1.\end{eqnarray} \tag{ 13 }$$

In particular, when the fundamental frequencies become commensurable, high-order or secondary spin–orbit resonances may happen. The condition of a certain k₁:k₂ spin–orbit resonance takes the form,

$$\begin{eqnarray}&&{k}_{1}{\dot{\sigma }}_{1}^{* }-{k}_{2}{\dot{\sigma }}_{2}^{* }=0\Rightarrow {k}_{1}{\dot{\sigma }}_{1}^{* }-{k}_{2}=0.\end{eqnarray} \tag{ 14 }$$

By solving Equation (14), it is possible to determine the nominal location of high-order and/or secondary resonances.

In Figure 3, the nominal locations of high-order and secondary spin–orbit resonances are shown in the phase space $(\theta ,\dot{\theta })$ for the spin–orbit problems with α = 0.35 and α = 0.65, and they are compared with the associated Poincaré sections. In both panels of Figure 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 those 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. By 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 Figure 4, we further present the distribution of high-order and secondary resonance in the $(\dot{\theta },\alpha )$ space and compare them to the dynamical map of the normalized second-derivative-based index ∣∣ΔD∣∣. The map of ∣∣ΔD∣∣ shown in the right-hand panel is the same as the one shown in Figure 2. In the left-hand 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 high-order 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 $\dot{\theta }=1$.

Zoom In Zoom Out Reset image size

From the right-hand panel of Figure 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∣∣.

4.2. Resonant Hamiltonian Model

From the previous subsection, it is known that high-order 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 reorganized as

$$\begin{eqnarray}&&{ \mathcal H }={{ \mathcal H }}_{0}\left({{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right)+{{ \mathcal H }}_{1}\left({\sigma }_{1}^{* },{\sigma }_{2}^{* },{{\rm{\Sigma }}}_{1}^{* },{{\rm{\Sigma }}}_{2}^{* }\right),\end{eqnarray} \tag{ 15 }$$

where the perturbation part ${{ \mathcal H }}_{1}$ is much smaller than the Kernel Hamiltonian ${{ \mathcal H }}_{0}$. Therefore, it is reasonable to take ${{ \mathcal 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 ${\sigma }_{1}^{* }$ and ${\sigma }_{2}^{* }$.

To study the k₁:k₂ resonance between ${\sigma }_{1}^{* }$ and ${\sigma }_{2}^{* }$, let us introduce the following transformation,

$$\begin{eqnarray}&&\begin{array}{l}{\gamma }_{1}={\sigma }_{1}^{* }-\displaystyle \frac{{k}_{2}}{{k}_{1}}{\sigma }_{2}^{* },\quad {{\rm{\Gamma }}}_{1}={{\rm{\Sigma }}}_{1}^{* },\\ {\gamma }_{2}={\sigma }_{2}^{* },\quad {{\rm{\Gamma }}}_{2}={{\rm{\Sigma }}}_{2}^{* }+\displaystyle \frac{{k}_{2}}{{k}_{1}}{{\rm{\Sigma }}}_{1}^{* },\end{array}\end{eqnarray} \tag{ 16 }$$

which is canonical with the generating function

$$\begin{eqnarray*}&&S={\sigma }_{1}^{* }{{\rm{\Gamma }}}_{1}+{\sigma }_{2}^{* }\left({{\rm{\Gamma }}}_{2}-\displaystyle \frac{{k}_{2}}{{k}_{1}}{{\rm{\Gamma }}}_{1}\right).\end{eqnarray*}$$

The resonant argument for the k₁:k₂ resonance is ${\gamma }_{1}={\sigma }_{1}^{* }-\tfrac{{k}_{2}}{{k}_{1}}{\sigma }_{2}^{* }$.

Under the new set of canonical variables given by Equation (16), the Hamiltonian can be further expressed as

$$\begin{eqnarray}&&{ \mathcal H }={{ \mathcal H }}_{0}\left({{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}\right)+{{ \mathcal H }}_{1}\left({\gamma }_{1},{\gamma }_{2},{{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}\right).\end{eqnarray} \tag{ 17 }$$

In particular, when the triaxial satellite is located inside a certain k₁:k₂ spin–orbit resonance, the corresponding resonant argument ${\gamma }_{1}={\sigma }_{1}^{* }-\tfrac{{k}_{2}}{{k}_{1}}{\sigma }_{2}^{* }$ becomes a long-periodic angle and the other one γ₂ is a short-periodic angle. This 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₁ times the period of γ₂, leading to the resonant Hamiltonian,

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal H }}^{* }=\displaystyle \frac{1}{2{k}_{1}\pi }{\displaystyle \int }_{0}^{2{k}_{1}\pi }{ \mathcal H }\left({\gamma }_{1},{\gamma }_{2},{{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}\right){\rm{d}}{\gamma }_{2}\\ =\,{{ \mathcal H }}_{0}\left({{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}\right)+\displaystyle \frac{1}{2{k}_{1}\pi }{\displaystyle \int }_{0}^{2{k}_{1}\pi }{{ \mathcal H }}_{1}\left({\gamma }_{1},{\gamma }_{2},{{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}\right){\rm{d}}{\gamma }_{2}.\end{array}\end{eqnarray} \tag{ 18 }$$

Under the resonant Hamiltonian model determined by ${{ \mathcal H }}^{* }$, the angle γ₂ becomes a cyclic variable, so that its conjugate momentum

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{2}={{\rm{\Sigma }}}_{2}^{* }+\displaystyle \frac{{k}_{2}}{{k}_{1}}{{\rm{\Sigma }}}_{1}^{* }\end{eqnarray*}$$

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 ${{ \mathcal H }}^{* }$ in the phase space (i.e., phase portraits). To compare phase portraits with Poincaré sections, we need to project the phase space $({\sigma }_{1}^{* },{{\rm{\Sigma }}}_{1}^{* })$ to the space $(\theta ,\dot{\theta })$ by considering the condition of Poincaré sections: $\mathrm{mod}\left({q}_{2},2\pi \right)=0$.

4.3. Results

Analytical and numerical results for high-order and secondary spin–orbit resonances in the phase space $(\theta ,\dot{\theta })$ are presented in Figure 5. Regarding analytical result, we provide the dynamical separatrix, which stands for the isolines of resonant Hamiltonian. The 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-hand column panels, dynamical separatrices of the synchronous primary resonance, high-order resonances, and secondary resonances are shown. In the right-hand column panels, the 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.

Zoom In Zoom Out Reset image size

From Figure 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:

$$\begin{eqnarray*}&&\left\{\begin{array}{l}\theta \left(t={t}_{0}\right)=\theta \left(t={t}_{0}+T\right),\\ \dot{\theta }\left(t={t}_{0}\right)=\dot{\theta }\left(t={t}_{0}+T\right).\end{array}\right.\end{eqnarray*}$$

Without loss of generality, we assume the angle at t = t₀ is equal to zero, i.e., θ₀ = 0. Thus, there are only two variables $({\dot{\theta }}_{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.

Figure 6 shows the analytical characteristic curves in the left-hand panel and comparisons between analytical and numerical curves in the right-hand panel. The shaded region stands for the space where the synchronous primary resonance happens. From Figure 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.

Zoom In Zoom Out Reset image size

Regarding the spin–orbit problem, Gkolias et al. (2016, 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 Figure 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 Figure 6).

Figure 7 shows all the analytical characteristic curves of the synchronous primary resonance, secondary resonances, and high-order resonances in the left-hand panel and makes a comparison between analytical curves to the numerical map of ∣∣ΔD∣∣ in the right-hand 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.

Zoom In Zoom Out Reset image size

In particular, Figure 7 provides a global view for the rotational dynamics of the secondary in the $(\dot{\theta },\alpha )$ 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.

The analytical approach is applied to two representative binary asteroid systems: (65803) Didymos and (4383) Suruga (Warner 2013; Pravec et al. 2016). (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 location of these two systems are marked in the right-hand panel of Figure 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⁻³ for (65803) Didymos and L_spin/L_orb = 9.9 × 10⁻⁴ 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 Figure 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 lines and the separatrix of the secondary 1:1 resonance is given in blue lines. From Figure 8, we can observe that (a) there is an excellent agreement between the analytical and numerical structures, and (b) the chaotic layer becomes wilder with a higher eccentricity.

Zoom In Zoom Out Reset image size

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 is 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 α smaller than ∼1.06. Thus, more observations are required to further determine the asphericity parameter α accurately to predict their real resonant states. However, Figure 7 can help to predict other possible resonant states when α is varied.

In this work, high-order and secondary resonances in the spin–orbit problem are investigated. Dynamical maps are numerically produced to provide a global picture for the dynamics. In particular, two types of dynamical indicators are introduced. The first is the fast Lyapunov indicator (FLI) and the second is the normalized second-derivative-based index ∣∣ΔD∣∣. The maps of FLI can distinguish the regular and chaotic regions, and the maps of ∣∣ΔD∣∣ have 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 those 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), a resonant Hamiltonian model is formulated. Thus, it is ready to study high-order and secondary resonance dynamics analytically. The 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 may be inside the secondary 1:1 spin–orbit resonance.

The main results are presented in Figure 7, where the distributions of synchronous primary resonance, high-order resonance and secondary resonances are presented and compared with dynamical map of ∣∣ΔD∣∣. This provides a global picture for us to predict possible resonant states when the asphericity parameter is changed.

The author wishes to thank an anonymous reviewer for helpful suggestions that improved the quality of this manuscript. 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).