© The Authors 2024

1 Introduction

Resonance lock with dissipative effects occurs in various natural systems, such as in the dynamics of dust (Dermott et al. 1994), the satellites of major planets (Goldreich & Peale 1966), and exoplanets (Ferraz-Mello et al. 2003). We study in this paper the motions inside the Jovian 2:1 mean motion resonance (abbreviated as J2:1) with the Yarkovsky effect (Peterson 1976; Farinella et al. 1998), one of the dissipative effects that can cause secular orbital evolution of asteroids. Outside the resonance, the Yarkovsky mobility of the asteroid is mainly the secular drift of the semi-major axis (Bottke et al. 2001). We focus on the Yarkovsky-driven secular evolution of eccentricities of asteroids inside the J2:1. The stimulation of eccentricity of resonant asteroids is closely related to the depletion of the asteroids inside the Kirkwood gaps, among which the J2:1 is known as the Hecuba gap, and also closely related to the formation of near-Earth asteroids, whose size-frequency (Morbidelli & Vokrouhlický 2003), orbit (Granvik et al. 2018), and spin-state (La Spina et al. 2004; Liu et al. 2023) distributions may reflect the selection effects of the Yarkovsky effect.

Being a gap in the man belt with the largest width, the dynamics of the J2:1 resonance has been studied in detail (Morbidelli 2002; Nesvorný et al. 2002). The overlap of secondary resonances was found to produce a zone of chaotic motions at low to moderate eccentricity (Henrard & Lemaitre 1987), while the overlap of secular resonances was found to produce a large chaotic zone at a large eccentricity (Morbidelli & Moons 1993). For a large J2:1 libration amplitude and moderate inclination, a bridge connecting the two chaotic regions was also found to provide a route to escape with initially small eccentricity (Henrard et al. 1995). Nevertheless, there remains two quasi-stable islands away from any low-order resonances (Moons et al. 1998), inside which slow diffusive chaos was detected using numerical methods (Ferraz-Mello et al. 1998b) and explained by the resonance with the great inequality (Ferraz-Mello et al. 1998a).

The Hecuba gap is not completely depleted of asteroids (Roig et al. 2002): 374 resonant J2:1 asteroids with diameters 1~20 km were confirmed by Chrenko et al. (2015), of which around 230 are long-lived asteroids (with dynamical lifetime >0.7 Gyr) that reside in the quasi-stable islands mentioned above and in their vicinity and around 140 are shorted-lived asteroids (with a dynamical lifetime less than 0.7 Gyr) that reside in the chaotic regions. Brož et al. (2005) showed that due to the Yarkovsky drift of the semi-major axis, prograde rotators could be continuously pushed into the region occupied by the short-lived J2:1 population in order to maintain a quasi-steady state, but they can rarely enter the quasi-stable islands. The origin of the long-lived asteroids, whether primordial or through capture from main belts, still remains unsolved (Chrenko et al. 2015).

The stability of stable resonant J2:1 asteroids with respect to the Yarkovsky effect was investigated by Brož & Vokrouhlický (2008) using N-body simulation incorporated with the Yarkovsky effect. The systematic drift of eccentricity inside the resonance – a common characteristic for resonance lock with dissipative effects (Gomes 1995) – was found; the dynamical lifetime as a function of asteroid size for prograde rotators to gain critical eccentricity was obtained; and retrograde rotators were shown to be more likely to escape. The pioneering work by Brož et al. on the systematic drifts of eccentricity of resonant asteroids was mainly carried out in a numerical way. In this paper, we carry out a perturbative treatment using the semi-analytical approach suitable for the J2:1. In addition to the interaction between the Yarkovsky effect with the J2:1 resonance, we further consider its interaction with the v₅ resonance inside the J2:1.

The contribution of this paper is twofold. The first part is a revisit of the resonance structures inside the J2:1 commensurability under the model of the planar elliptic restricted three-body problem (ERTBP). The first part has three components: the Schubart averaging of the planar ERTBP, Morbidelli’s elimination of the J2:1 harmonic (to be more specific, semi-analytical computation of the Arnold action variable of J2:1 resonance and the coordinate transform by the Henrard-Lemaitre formula), and the partial normalization of the secular Hamiltonian for the v₅ resonance to the fourth order. For the first component, we complement the method developed in Moons (1994) by computing the transformation of variables between the complete system and the averaged system. The third component can be considered as a slight extension of the treatment of v₅ resonance in Morbidelli & Moons (1993), where the integrable approximation of v₅ resonance is obtained by a first-order perturbation theory.

The second contribution, which to our knowledge is original, is a perturbative treatment of the Yarkovsky-driven systematic drift in the action space of the J2:1 resonance. Under the model of planar ERTBP, we first identify two patterns of drift in the action space through numerical simulations: the first one (Type I) – already pointed out by Brož & Vokrouhlický (2008) – “walks” in the action space almost parallel to the curve of the stable resonance center of the J2:1 resonance; the second one (Type II) “walks” almost parallel to the curve of the stable resonance center of the v₅ resonance. Our approach to analyze these two types of drift is close to that of Gomes (1995): the suitable form of the Yarkovsky effect averaged over the mean motions and the J2:1 resonance cycle is applied to the integrable approximations of the J2:1 resonance and the v₅ resonance in order to obtain quasi-Hamiltonian systems from which a vector field of drift rates in the action space is derived.

For the Type I drift, the eccentricity increases (decreases) along the vector field for prograde (retrograde) asteroids. For prograde asteroids, we propose a fast algorithm to compute the Yarkovsky-driven excitation of the orbit eccentricity. Moreover, the Type I drift with small eccentricities and small libration amplitudes of the J2:1 resonance can be modeled by a harmonic oscillator with a slowly varying parameter. An analytical treatment of this simplified model using the classical adiabatic invariant theory was carried out and shows that for prograde (retrograde) asteroids, the J2:1 libration amplitude decreases (increases). We also find that the Yarkovsky-driven Type I drift for retrograde asteroids may serve as the mechanism to drive asteroids to the bridge identified by Henrard et al. (1995) that transfers them from the low-eccentricity chaotic regions caused by overlap of secondary resonances to the high-eccentricity chaotic regions caused by overlap of secular resonances. The N-body simulation results incorporated with the Yarkovsky effect confirm the existence of this evolutionary route for retrograde asteroids to escape the J2:1 resonance.

This paper is organized as follows. Section 2 presents the model and the method. Section 3 presents the results of the first part of the study on the resonance dynamics in the J2:1 resonance. Section 4 presents the results of the second part of the study on the Yarkovsky-driven drift in the action space. Section 5 discusses the Yarkovsky-driven escape of retrograde asteroids from the J2:1, and Sect. 6 concludes the study.

Fig. 1 Schematic procedures of the perturbative treatments for the Yarkovsky-driven drifts in the J2:1 commensurability under the model of planar ERTBP.

2 Model and method

The method adopted in this paper is schematically presented in Fig. 1. It consists of two chains corresponding to the two parts of studies mentioned in the introduction: (1) semi-analytical investigations on the dynamical structures in the J2:1 commen-surability in the canonical formalism. That is, a chain of three canonical transformations are applied to obtain a normal form of the Hamiltonian which can approximate well the secular dynamics inside the J2:1 commensurability; (2) averaging of the Yarkovsky force based on the aforementioned canonical transformations. The averaged Yarkovsky forces can be inserted into the canonical equations of the Hamiltonian system at different levels.

The first part starts with the Hamiltonian function H_pERTBP of the planar elliptic restricted three-body problem (ERTBP). We include in the model a single-frequency precession of the Jupiter’s mean longitude of perihelion. Adopting a suitable canonical chart (P, Q) for the J2:1 resonance, we consult the Schubart averaging to average the Hamiltonian over the mean motions to obtain an averaged Hamiltonian ${\overline{H}}_{\text{Schubart }}(\overline{P},\overline{Q})$, which is bridged with H_pERTBP (P, Q) by the canonical transformation ℒ₁. Then, a 1-degree-of-freedom (1-DOF) integrable approximation ${\overline{H}}_a$ which describes the J2:1 resonance is extracted from ${\overline{H}}_{\text{Schubart}}$ by averaging. Next, we introduce the Arnold action variable of ${\overline{H}}_a$, and develop the Fourier expansion of the Hamiltonian in terms of the new chart (I, θ), denoted by ${\overline{H}}_{\text{Fourier}}$, which is bridged with ${\overline{H}}_{\text{Schubart}}$ by ℒ₂. At last, we apply a Lie transform ℒ₃ which transforms the Hamiltonian into a partial normal form ${\overline{H}}_{\text{Normal }}^{\ast }$ (J, ψ), which describes the v₅ resonance inside the J2:1 commensurability.

The second part starts with the linear analytical model of diurnal Yarkovsky effect by Vokrouhlický (1998). The Yarkovsky force is originally expressed in terms of position and velocity. To insert the Yarkovsky force into the canonical equation of H_pERTBP, we first rewrite it in terms of the chart (P, Q) as F_Y|(P,q). Then, suitable averaging of the Yarkovsky force is carried out using the above chain of canonical transformations.

2.1 Hamiltonian of planar elliptic restricted three-body problem with the second body’s orbit in precession

We adopted the planar ERTBP to model the three-body system Sun(m₀)-asteroid(m₁)-Jupiter(m₂). Denote (a_i, e_i, M_i,${\widetilde{ɡ}}_i$, λ_i) as the semi-major axis, the eccentricity, the mean anomaly, the mean longitude of perihelion, and the mean longitude of body m_i, with i = 1,2. The Hamiltonian function of the planar ERTBP is $$H=-\frac{{\mu }_1}{2a_1}+\left(n_2+g_5\right){\text{Λ}}_2+g_5{\text{Γ}}_2+G{m}_2\left(\frac{r_1\cdot r_2}{r_2^3}-\frac{1}{\left|r_1-r_2\right|}\right),$$(1)

where µ₁ = G (m₀ + m₁) = Gm₀ with G the gravitational constant; r_i is the position vector from m₀ to m_i, and r_i = |r_i|. We assume that m₂ moves on a Keplerian orbit around the Sun, which undergoes a uniform mean motion Ṁ₂ = n₂ = $\sqrt{G\left(m_0+m_2\right)/a_2^3}$, and a precession of mean longitude of perihelion ${\dot{\widetilde{ɡ}}}_2={ɡ}_5$. To make the Hamiltonian autonomous, two extra degrees of freedom (λ₂, Λ₂) and $\left({\widetilde{ɡ}}_2,{\text{Γ}}_2\right)$ are introduced, which together with the two degrees of freedom of the asteroid’s planar motion make the Hamiltonian Eq. (1) a 4-DOF system.

To describe the motion of an asteroid trapped in the p + q : p mean motion resonance with Jupiter, we adopt the canonical variables $$\{\begin{array}{ll}{P}_1=-\frac{q}{p}{L}_1,\hfill & Q_1=\frac{p+q}{q}{\lambda }_2-\frac{p}{q}{\lambda }_1-{\widetilde{ɡ}}_1,\hfill \\ P_2=\left(\frac{q}{p}+1\right)L_1-G_1,\hfill & Q_2={\widetilde{ɡ}}_2-{\widetilde{ɡ}}_1,\hfill \\ P_3=\frac{p+q}{p}{L}_1+{\text{Λ}}_2,\hfill & Q_3={\lambda }_2,\hfill \\ P_4={\text{Γ}}_2-P_2,\hfill & Q_4={\widetilde{ɡ}}_2,\hfill \end{array}$$(2)

where $L_1=\sqrt{{\mu }_1{a}_1}$, and $G_1=L_1\sqrt{1-e_1^2}$. In terms of the chart Eq. (2), the Hamiltonian Eq. (1) is rewritten as $$\begin{array}{c}H(P,Q)=-{\left(\frac{q}{p}\right)}^2\frac{{\mu }_1^2}{2P_1^2}+n_2\frac{p+q}{q}{P}_1+n_2{P}_3+{ɡ}_5\left(P_2+P_4\right)\\ +G{m}_2\left(\frac{r_1\cdot r_2}{r_2^3}-\frac{1}{\left|r_1-r_2\right|}\right).\end{array}$$(3)

In this paper, we focus on 2:1 commensurability such that p = q = 1. The system parameters of the Hamiltonian are taken as m₂/m₀ = 10⁻³, e₂ = 0.04419, and ɡ₅ = 0¹. The normalized units are adopted: the unit of mass is [M] = m₀ + m₂ = 1.001 × 1.98910³⁰ kg, the unit of Length [L] = a₂ = 5.2 AU, and the unit of time $[T]=\sqrt{{[L]}^3/G[M]}\approx 1.88632\text{\hspace{0.17em}yr}$. Under this setting we have G = 1[L]³[M]⁻¹[T]⁻² and n₂ = 1[T]⁻¹.

2.2 Schubart averaging of the planar elliptic restricted three-body problem

The first step to reduce the number of degrees of freedom of the system is to average the Hamiltonian over the mean motions of the system which have the shortest timescales and thus are well separated from other degrees of freedom (e.g., Beaugé & Roig 2001).

Using the canonical variables Eq. (2) and the normalized units above, Q₃ which represent the mean motions of the system has a period 2π. By the inversion of Eq. (2), the angle Q₃ appears in H(P,Q) in the form of Q₃/p. So, by the averaging theorem, we may average the Hamiltonian function over Q₃ to obtain the averaged system: $$\overline{H}(\overline{P},\overline{Q})=\frac{1}{2p\pi }{{\displaystyle {\displaystyle \underset{0}{\overset{2p\pi }{\int }}H\left(P,Q\right)\text{d}{Q}_3.}}}^{\text{}}$$(4)

It is important to note that according to the d’Alembert principle, Q₃ must appear in combinations with Q₄. So, the averaged system $\overline{H}$ is independent of both ${\tilde{Q}}_3,{\tilde{Q}}_4,$ which makes ${\overline{P}}_3,{\overline{P}}_4$ two integrals of motion. $\overline{H}$ is thus a 2-DOF system.

To compute the integral Eq. (4) analytically, it is customary to expand the Hamiltonian into a Fourier series in Q_i; however, this is only valid at small eccentricity due to the convergence problem of the analytical expansion method. To avoid this problem, we consulted the numerical averaging method, also known as the Schubart averaging (Moons 1994). The idea is to replace the integration in Eq. (4) by the Riemannian sum, namely² $$\begin{array}{c}\overline{H}(\overline{P},\overline{Q})=\frac{1}{N}{\displaystyle \sum _{n=1}^{N}H\left(\begin{array}{c}{Q}_1={\overline{Q}}_1,Q_2={\overline{Q}}_2,Q_3=\frac{2np\pi }{N},\\ P_1={\overline{P}}_1,P_2={\overline{P}}_2,P_3={\overline{P}}_3\end{array}\right)}\\ +{\epsilon }_1(N;\overline{P},\overline{Q}).\end{array}$$(5)

The canonical equation of $\overline{H}(\overline{P},\overline{Q})$, which cannot be derived directly since $\overline{H}(\overline{P},\overline{Q})$ is not known explicitly, is obtained by conducting Schubart averaging of the canonical equation of H (P, Q). Details for this can be found in Appendix A.1. We note that though the formulae presented in Appendix A.1 are singular at e₁ = 0 because of the usage of Kepler elements, it is sufficient for our studies in this paper. For nonsingular formulae, one may consult Moons (1994) and references there.

A remark is that the result of direct averaging obtained in Eq. (4) is identical with the first-order result of standard perturbation method (such as Lie transform). In Appendix A.2, we present the Lie mapping $${ℒ}_1:(\overline{P},\overline{Q})\to (P,Q),$$(6)

which bridges H(P, Q) and $\overline{H}(\overline{P},\overline{Q})$ (accurate to the first order). To our knowledge, the computation of ℒ₁ using Schubart averaging has not been shown in previous researches (though the necessity of this transformation has been mentioned by Moons 1994).

2.3 Morbidelli’s elimination of J2:1 harmonic

The consequence of the Schubart averaging result $\overline{H}(\overline{P},\overline{Q})$ not being known explicitly is that extra efforts are required for the application of the standard perturbation method. We follow the method developed by Morbidelli (1993) in the following.

Suppose $\overline{H}(\overline{P},\overline{Q})$ is expanded into a Fourier series of ${\overline{Q}}_1$ and ${\overline{Q}}_2$, the J2:1 harmonics, namely the trigonometric terms only containing the phase ${\overline{Q}}_1$, would have the largest strength³ for resonant asteroids in the J2:1. So, by simply removing all harmonics involving ${\overline{Q}}_2$, we can extract from $\overline{H}$ an integrable approximation ${\overline{H}}_a$ for the J2:1 resonance as $${\overline{H}}_a\left({\overline{Q}}_1,{\overline{P}}_1;{\overline{P}}_2\right)=\frac{1}{2\pi }{{\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}\overline{H}(\overline{P},\overline{Q})\text{d}{\overline{Q}}_2.}}}^{\text{}}$$(7)

The integrable approximation ${\overline{H}}_a$ is 1-DOF with ${\overline{P}}_2$ as a parameter and could describe well the motions inside the J2:1 within a time interval comparable to the period of the J2:1 resonance. We note that the averaging in Eq. (7) removes a remainder term $R(\overline{P},\overline{Q}):=\overline{H}-{\overline{H}}_a$.

Because $\overline{H}$ is not known explicitly, the integration in Eq. (7) is also replaced by the Riemann summation $$\begin{array}{c}{\overline{H}}_a\left({\overline{Q}}_1,{\overline{P}}_1;{\overline{P}}_2\right)\\ =\frac{1}{N}{\displaystyle \sum _{n=1}^N\overline{H}(\begin{array}{c}{Q}_1={\overline{Q}}_1,{\overline{Q}}_2=\frac{2n\pi }{N},\\ P_1={\overline{P}}_1,P_2={\overline{P}}_2\end{array})+{\epsilon }_2(N;\overline{P},\overline{Q})}.\end{array}$$(8)

Next, we use Morbidelli’s elimination of harmonic to eliminate the J2:1 harmonic in the system $\overline{H}(\overline{P},\overline{Q})$, that is, to apply a change of variable to the system $\overline{H}$ such that ${\overline{H}}_a$ is independent of phase variables. The formal procedure of the method consists of two steps: (1) the numerical computation of the action-angle (I₁,θ₁) variable of the integrable approximation ${\overline{H}}_a$, which corresponds to a canonical transformation ${ℒ}_2^{\prime }:\left(I_1,{\theta }_1\right)\to \left({\overline{P}}_1,{\overline{Q}}_1\right)$, and (2) the extension of ${ℒ}_2^{\prime }$ to a canonical transformation $${ℒ}_2:(I,\theta )\to (\overline{P},\overline{Q}),$$(9)

which applies to the full system $\overline{H}$. Under ℒ₂, the J2:1 harmonic is eliminated in the sense that ${\overline{H}}_a\circ {ℒ}_2(I,\theta )$ becomes independent of θ₁, θ₂, namely ${\overline{H}}_a\circ {ℒ}_2(I,\theta )={\overline{H}}_a(I)$. Moreover, the Fourier expansion of the Hamiltonian $\overline{H}$ under ℒ₂ becomes $$\begin{array}{c}\overline{H}(I,\theta )={\overline{H}}_a\left(I_1,I_2\right)\\ +{\displaystyle \underset{\begin{array}{c}{k}_1\ge 0,k_2\in Z/\left\{0\right\}\\ \gcd \left(k_1,k_2\right)=1\end{array}}{{\sum }^{\text{}}}}{\displaystyle \sum _{l=1}^{\infty }{c}_{k_1,k_2}^l\left(I_1,I_2\right)\cos \left(l\left(k_1{\theta }_1+k_2{\theta }_2\right)\right).}\end{array}$$(10)

The above formal description of Morbidelli’s method to eliminate a target harmonic in the Hamiltonian may appear as standard at first. However, the real difficulty in practice is that $\overline{H}(\overline{P},\overline{Q})$ and ${\overline{H}}_a(\overline{P},\overline{Q})$ are not known explicitly. The method developed by Morbidelli (1993) allows one to obtain the values of the functions ℒ₂, ${\overline{H}}_a$, and $c_{k_1,k_2}^l$ on a tabulated grid in the action space (I₁, I₂) (see the grid 𝒢₃defined in Eq. (B.10)) Then, the value of $\overline{H}$ (I, θ) as well as its partial derivatives can be obtained by interpolation and numerical differentiation on any point inside the grid. Details on the application of Morbidelli’s method to compute ℒ₂ and $\overline{H}$ (I, θ) are tedious, and saved in Appendices B.2 and B.3.

2.4 Partial normalization of secular Hamiltonian for v₅ resonance

With the Fourier series Eq. (10) at hand, following the standard perturbation theory (Ferraz-Mello 2007), we could use ${\overline{H}}_a$ as the kernel function to compute the normal form of the Hamiltonian. However, though ${\overline{H}}_a$ depends both on I₁, I₂, the system is degenerate in the sense that ${\dot{\theta }}_1=\left|\partial {\overline{H}}_a/\partial I_1\right|\gg {\dot{\theta }}_2={\left|\partial {\overline{H}}_a/\partial I_2\right|}^{}$⁴. In fact, for flows trapped in the v₅ resonance, ${\overline{H}}_a$ is not a valid integrable approximation for a full normalization. So, in the following, we use ${\overline{H}}_a$ as the kernel function to partially normalize the system, that is, a Lie transform of order n $${ℒ}_{3,n}:\left(J_1,{\psi }_1,J_2,{\psi }_2\right)\to \left(I_1,{\theta }_1,I_2,{\theta }_2\right),$$(11)

is applied to remove all harmonics involving θ₁ in $\overline{H}$ (I, θ) so that the transformed system, referred to as the partial normal form ${\overline{H}}_n^{\ast }$, takes the form $$\begin{array}{c}{\overline{H}}_n^{\ast }\left(J_2,{\psi }_2;J_1\right)=\overline{H}\circ {ℒ}_3(J,\psi )\\ =K_0(J)+K_1\left(J,{\psi }_2\right)+\dots +K_n\left(J,{\psi }_2\right)+R_{n+1}\left(J,{\psi }_1,{\psi }_2\right),\end{array}$$(12)

where $$\{\begin{array}{l}{K}_0\left(J_1,J_2\right)={\overline{H}}_a\left(J_1,J_2\right)\hfill \\ K_1\left(J_1,J_2,{\psi }_2\right)={\displaystyle \sum _{l=1}^{\infty }{c}_{0,1}^l\left(J_1,J_2\right)\cos \left(l{\psi }_2\right),}\hfill \\ \dots \dots \hfill \end{array}$$(13)

Ignoring the remainder R_n+1 (J, ψ), the partial normal form ${\overline{H}}_n^{\ast }$ is independent of the phase ψ₁ and contains only the phase angle ${\psi }_2\approx {\widetilde{ɡ}}_1-{\widetilde{ɡ}}_2$, that is, the phase angle of the v₅ resonance. So, ${\overline{H}}_n^{\ast }$ is in fact an integrable approximation of the v₅ resonance.

The above normalization procedure is again standard while the real difficulty in practice is that the Hamiltonian function $\overline{H}$ (I, θ) is not known explicitly, but it can be evaluated by interpolation on any point inside the tabulated grid. For details on the formal and the practical computations of the Lie mapping ℒ_3,n and the partial normal form ${\overline{H}}_n^{\ast }$, one may refer to Appendix B.4.

We end by remarking that the semi-analytical normalization algorithm is restricted to low orders (such as order three or four) since the higher order normal form requires higher order derivatives of the original Hamiltonian, which, generally, cannot be provided by standard interpolation methods.

2.5 Averaging of the diurnal Yarkovsky effect

We recall from Sect. 2.1 that the planar configuration is assumed for orbital motion in this paper. We further assume here that the spin-axis of the asteroid is normal to its orbital plane. In this case, the seasonal Yarkovsky effect (Vokrouhlický & Farinella 1998) vanishes and the diurnal Yarkovsky effect only has in-plane components. Consequently, the planar configuration persists throughout the evolution.

In our study, we adopt the classical analytical model of the diurnal Yarkovsky effect developed in Vokrouhlický (1998). In Appendix D, the formulae of the diurnal Yarkovsky force adopted in our study are presented. In particular, Eq. (D.6) presents the Yarkovsky force in terms of the chart (P, Q) in Eq. (2), denoted as F_Y|_(P,Q). Therefore, the task is to derive the averaged Yarkovsky force in terms of new charts generated by the canonical transformations in the previous sections.

Before proceeding, it is important to note that the correctness of the canonical equation of H (P, Q) with g₅ = 0 perturbed by F_Y|_(P,Q) can be checked using the augmented version of the SWIFT package which includes the classical analytical model of the Yarkovsky effect (Levison & Duncan 1994; Brož et al. 2011). In this way, we use the integration results of SWIFT package as the reference to verify the results obtained by the normal form of the Hamiltonian and the corresponding averaged Yarkovsky force.

2.5.1 Case of $\overline{H}(\overline{P},\overline{Q})$

To derive the averaged Yarkovsky force for $\overline{H}(\overline{\mathit{P}},\overline{\mathit{Q}})$, it is natural to average F_Y |_(P,Q) over one period of the phase Q₃, namely $${{\overline{F}}_Y|}_{(\overline{P},\overline{Q})}=\frac{1}{2p\pi }{\displaystyle {\int }_0^{2p\pi }{F}_Y}(P,Q)\text{d}{Q}_3.$$(14)

We recall that the averaged system $\overline{H}(\overline{\mathit{P}},\overline{\mathit{Q}})$ is 2-DOF. So the averaged Yarkovsky force ${{\overline{\mathit{F}}}_Y|}_{(\overline{\mathit{P}},\overline{\mathit{Q}})}$ includes only the differential rate for ${\overline{P}}_1,{\overline{P}}_2,{\overline{Q}}_1,{\overline{Q}}_2$, while formulae for ${\overline{P}}_3,{\overline{P}}_4,{\overline{Q}}_3,{\overline{Q}}_4$ are not needed. We note that the integration of Eq. (14) is also replaced by a Riemann summation, as in Eq. (5).

2.5.2 Case of ${\overline{H}}_a$

To insert the Yarkovsky force into ${\overline{H}}_a$, one conducts an averaging of ${{\overline{\mathit{F}}}_Y|}_{(\overline{\mathit{P}},\overline{\mathit{Q}})}$ over ${\overline{Q}}_2$ similar to Eq. (7). However, we note that the diurnal Yarkovsky effect is in nature local in the sense that it depends only on the radial distance of the asteroid from the Sun and not on the relative geometry with respect to the perturbing body Jupiter. So ${{\overline{\mathit{F}}}_Y|}_{(\overline{\mathit{P}},\overline{\mathit{Q}})}$ which is independent of the phase Q₃, should also be independent of all phase angles, that is, ${{\overline{\mathit{F}}}_Y|}_{(\overline{\mathit{P}},\overline{\mathit{Q}})}={{\overline{\mathit{F}}}_Y|}_{(\overline{\mathit{P}},\overline{\mathit{Q}})}(\overline{\mathit{P}})$. Thus, ${\overline{\mathit{F}}}_Y$ can be directly applied to the canonical equation of ${\overline{H}}_a$.

2.5.3 Case of $\overline{H}(\mathit{I},\mathit{\theta })$

The next is to find trie Yarkovsky force in terms of the chart (I, θ), denoted as ${{\overline{\mathit{F}}}_Y|}_{(I,\theta )}$. Using the canonical transformation ${ℒ}_2^{-1}:(\overline{\mathit{P}},\overline{\mathit{Q}})\to (\mathit{I},\mathit{\theta })$, it is natural to write $${{\overline{F}}_Y|}_{(I,\theta )}={\frac{\partial (I,\theta )}{\partial (\overline{P},\overline{Q})}{\overline{F}}_Y|}_{\overline{P},\overline{Q})}(\overline{P}).$$(15)

Because $I_2={\overline{P}}_2$ under ℒ₂ (see Eq. (B.6)), we have ${{\overline{\mathit{F}}}_Y|}_{I_2}={{\overline{\mathit{F}}}_Y|}_{{\overline{P}}_2}$. However, for the other three variables I₁,θ₁,θ₂, it is inconvenient to directly use Eq. (15) because ${ℒ}_2^{-1}$ is not known explicitly, but it requires numerical integration of the system ${\overline{H}}_a$ (see Eq. (B.8)). An alternative approach to compute ${{\overline{F}}_Y|}_{I_1}$ is as follows. We recall that the action variable I₁ of ${\overline{H}}_a$ corresponds to the oriented area of the phase flow of ${\overline{H}}_a$, and the phase flow of ${\overline{H}}_a$ is uniquely determined by its energy ${\overline{H}}_a=h$ since it is 1-DOF. So we have $I_1=I_1\left(h,{\overline{P}}_2\right)$, which yields $$\begin{array}{l}{{\overline{F}}_Y|}_{I_1}=\frac{\partial I_1\left(h,{\overline{P}}_2\right)}{\partial h}\left({\frac{\partial {\overline{H}}_a}{\partial {\overline{P}}_1}{\overline{F}}_Y|}_{{\overline{P}}_1}+\frac{\partial {\overline{H}}_a}{\partial {\overline{P}}_2}{\overline{F}}_{Y\mid {\overline{P}}_2}+{\frac{\partial {\overline{H}}_a}{\partial {\overline{Q}}_1}{\overline{F}}_Y|}_{{\overline{Q}}_1}\right)+\\ {\frac{\partial I_1\left(h,{\overline{P}}_2\right)}{\partial {\overline{P}}_2}{\overline{F}}_Y|}_{{\overline{P}}_2}.\end{array}$$(16)

The unknowns in Eq. (16), the two partial derivative $\partial I_1/\partial h,\partial I_1/\partial {\overline{P}}_2$, can be obtained by interpolation and numerical differentiation of the tabulated function I₁ (h, I₂) (see the end of Appendix B.2). It seems that there are no convenient ways to compute the differential rates of θ₁ , θ₂ except using Eq. (15), which requires too much computation. Considering that ${{\overline{\mathit{F}}}_Y|}_{\theta }$ is usually negligible comparing to the fundamental frequencies of θ₁ , θ₂, we simply ignore these two terms in the following of this paper, that is, we simply set ${\overline{\mathit{F}}}_{{Y|}_{{\theta }_1}}={{\overline{\mathit{F}}}_Y|}_{{\theta }_2}=0$. To conclude, we have $${{\overline{F}}_Y|}_I=\left(\begin{array}{c}{{\overline{F}}_Y|}_{I_1}\\ {{\overline{F}}_Y|}_{I_2}\end{array}\right)={G(\overline{P},\overline{Q}){\overline{F}}_Y|}_{(\overline{P},\overline{Q})}={\left(\begin{array}{cccc}\frac{\partial I_1}{\partial h}\frac{\partial {\overline{H}}_a}{\partial {\overline{P}}_1}& \frac{\partial I_1}{\partial h}\frac{\partial {\overline{H}}_a}{\partial {\overline{P}}_2}+\frac{\partial I_1}{\partial {\overline{P}}_2}& \frac{\partial I_1}{\partial h}\frac{\partial {\overline{H}}_a}{\partial {\overline{Q}}_1}& 0\\ 0& 1& 0& 0\end{array}\right){\overline{F}}_Y|}_{(\overline{P},\overline{Q})}.$$(17)

A serious inconvenience of Eq. (17) is that ${{\overline{\mathit{F}}}_Y|}_I$ is computed in terms of the chart $(\overline{\mathit{P}},\overline{\mathit{Q}})$ instead of the chart (I, θ). So to use them in the numerical integration of $\overline{H}(\mathit{I},\mathit{\theta })$ at every step, we need to first conduct a transformation ${ℒ}_2:(\mathit{I},\mathit{\theta })\to (\overline{\mathit{P}},\overline{\mathit{Q}})$. This requires heavy computations of multi-dimensional interpolations. One way to get around this is to conduct a further averaging over θ₁ as shown in Eq. (19) below.

2.5.4 Case of ${\overline{H}}^{*}(\mathit{J},\mathit{\psi })$

To incorporate the Yarkovsky force into ${\overline{H}}^{*}(\mathit{J},\mathit{\psi })$, we should conduct the following averaging: $$\begin{array}{l}{{\overline{F}}_Y|}_J(J,\psi )={\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\overline{F}}_Y}|}_I°{ℒ}_3(J,\psi )d{\psi }_1\\ ={\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }G}(\overline{P},\overline{Q}){\overline{F}}_Y|}_{(\overline{P},\overline{Q})}°{ℒ}_2°{ℒ}_3(J,\psi )d{\psi }_1,\end{array}$$(18)

where in the second equality we use Eqs. (9) and (17) so that the right-hand side is in terms of the chart $(\overline{\mathit{P}},\overline{\mathit{Q}})$ and can be directly calculated. Obviously, the averaging over ψ₁ implies that the result ${{\overline{\mathit{F}}}_Y|}_J(\mathit{J},\mathit{\psi })$ should only depend on ψ₂. The computation of Eq. (18) at every step of the numerical integration of ${\overline{H}}^{*}$ requires too much computation time (because of the computation of ℒ₂ as mentioned earlier). First observe that ℒ₃ (J, ψ) is a near-identity transformation. So we adopted the following simplification: $${{{\overline{F}}_Y|}_J(J,\psi )\approx \frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }G}{\overline{F}}_Y|}_{(\overline{P},\overline{Q})}°{ℒ}_2(I=J,\theta =\psi )d{\theta }_1.$$(19)

The result ${{\overline{\mathit{F}}}_Y|}_{\mathit{J}}(\mathit{J},\mathit{\psi })$ in Eq. (19) should not depend on ψ₂ either, that is, ${{\overline{\mathit{F}}}_Y|}_{\mathit{J}}={{\overline{\mathit{F}}}_Y|}_{\mathit{J}}(\mathit{J})$ (because the mapping result $\left({\overline{P}}_1,{\overline{P}}_2\right)$ given by ℒ₂ (J, ψ) does not depend on the value of ψ₂, and ${\overline{\mathit{F}}}_Y={\overline{\mathit{F}}}_Y\left({\overline{P}}_1,{\overline{P}}_2\right)$ as mentioned earlier). This prompts us to compute ${{\overline{\mathit{F}}}_Y|}_{\mathit{J}}$ on a grid of the action space (I₁, I₂). In fact, the computation of Eq. (19) becomes straightforward for each node of the grid 𝒢₃, where the mapping ℒ₂ is known (for this see Appendix B.3). In this way, we obtain a tabulated function ${\left({{\overline{\mathit{F}}}_Y|}_{\mathit{J}}\right)}_{i,j}$ on the grid 𝒢₃, and ${{\overline{\mathit{F}}}_Y|}_{\mathit{J}}(\mathit{J})$ can be evaluated by interpolation for any point inside the grid. For ℒ₃ (J, ψ) is near-identity, we take the simplification that ${{\overline{\mathit{F}}}_Y|}_{\mathit{I}}(\mathit{I})={{\overline{\mathit{F}}}_Y|}_{\mathit{J}}(\mathit{I})$ so that it can be applied to the system $\overline{H}(I,\theta )$ as well.

3 Dynamics in the J2:1 commensurability

In this section we present the dynamical structures inside the J2:1 commensurability, mainly the location, stability, width of the J2:1 and the ν₅ resonances. We show first how these structures are derived from the normal form of the Hamiltonian, and then how these structures depend on the normalization order. We also show by surfaces of section and the numerical integration of example orbits that these structures presented here reflect properly (both qualitatively and quantitatively) the dynamics of the original model of planar ERTBP, except for flows close to the J2:1 libration boundaries.

3.1 Resonance structures inside the J2:1 commensurability

3.1.1 Resonance structure of the integrable approximation ${\overline{H}}_a$ for the J2:1 libration

From Sect. 2 it is known that ${\overline{H}}_a\left({\overline{P}}_1,{\overline{Q}}_1;{\overline{P}}_2\right)$ is a 1-DOF Hamil-tonian system with a parameter ${\overline{P}}_2$. Analysis of the dynamics of this 1-DOF system is straightforward, and well developed in literatures. So, we save the technical details in Appendix B.1, and directly present the result here.

Figure 2 summarizes the resonance structures of ${\overline{H}}_a\left({\overline{P}}_1,{\overline{Q}}_1;{\overline{P}}_2\right)$. We note that instead of the $\left({\overline{P}}_1,{\overline{P}}_2\right)$ plane, the structures are presented on the ${\overline{a}}_1-{\overline{e}}_1$ plane using the relation Eq. (2), which is more intuitive. In Fig. 2, the red (blue) curve marks the locations of the equilibriums of ${\overline{H}}_a\left({\overline{P}}_1,{\overline{Q}}_1;{\overline{P}}_2\right)$ at ${\overline{Q}}_1=0\left({\overline{Q}}_1=\pi \right)$. The equilibrium at ${\overline{Q}}_1=0$ is always a center. The yellow and green curves mark the boundaries of the libration domain at the section ${\overline{Q}}_1=0$, that is, they consist of the intersection points of the largest libration orbits with the section ${\overline{Q}}_1=0$ (see Appendix B.1). Then, the libration domain for fixed ${\overline{P}}_2$ is uniquely characterized by the segment of the level curve of ${\overline{P}}_2$ (presented as horizontal gray dashed curves in Fig. 2) in between the libration boundaries. Any point on the segment corresponds uniquely to a libration flow. Also, we note that the shrink of the libration width with ${\overline{e}}_1>0.57$ is due to our numerical setting to avoid precision loss of the Schubart averaging, which suffers from the collision singularity at |r₁ – r₂| = 0 (see Appendix B.1 for details).

The situation of the equilibrium at ${\overline{Q}}_1=\pi$ is somewhat more complicated: (1) The equilibrium is absent at small eccentricity (e.g., see frame a of Fig. B.2); (2) The equilibrium at ${\overline{Q}}_1=\pi$ (if exist) is a saddle for small and moderate eccentricity (≤ 0.5), and becomes a center as eccentricity further increases (Moons & Morbidelli 1993); (3) The un-smoothness of the saddle curve at moderate eccentricity is due to the precision loss of the Schubart averaging, which originates from the collision singularity at |r₁ – r₂| = 0. The libration boundaries of the center at ${\overline{Q}}_1=\pi$ for large eccentricity are not presented.

From the above, each coordinate $\left({\overline{P}}_1,{\overline{P}}_2\right)$ in the libration domain in Fig. 2 uniquely characterizes a libration flow of ${\overline{H}}_a$ with initial condition $\left({\overline{P}}_1,{\overline{P}}_2,{\overline{Q}}_1=0\right)$. The inverse of the canonical transformation ℒ₂ in Eq. (9) maps $\left({\overline{P}}_1,{\overline{P}}_2,{\overline{Q}}_1=0,{\overline{Q}}_2\right)$ to the action variable (I, θ) where I1 corresponds to the area encircled by the libration flow (for details in see Appendix B.2). In Fig. 2, the level curves of I₁ are presented as vertical black solid curves. That is, the oriented area of the libration flow on the level curves of I1 remains constant. Particularly, the red curve for the resonance center is the level curve I₁ = 0.

We recall that ${\overline{H}}_a(I)={\overline{H}}_a(\overline{P},\overline{Q})°{ℒ}_2(I,\theta )$ where the phase angle θ₁ represents the J2:1 libration angle, and θ₂ represents the ν₅ resonance angle (see Eq. (2) and Eq. (B.8)). Ignoring the perturbations from all harmonics in Eq. (10), we obtain the two unperturbed frequencies ${\omega }_i=\partial {\overline{H}}_a/\partial I_i$ with i = 1,2. It is known that the overlap of the secondary resonances of the form ω₁ – kω₂ produces the chaotic motions at small eccentricity e₁ ≤ 0.2 in the J2:1 commensurability (Henrard & Lemaitre 1987). In Fig. 2, the colored curves with ${\overline{e}}_1<0.2$ present the locations of the secondary resonances with k = 2,…, 5.

Fig. 2 Resonance structures (equilibriums and libration width) of the integrable approximation ${\overline{H}}_a$ for the J2:1 resonance.

3.1.2 Resonance structure of the integrable approximation ${\overline{H}}^{*}$ for the ν₅ resonance

Similar to ${\overline{H}}_a\left({\overline{P}}_1,{\overline{Q}}_1;{\overline{P}}_2\right)$, the partial normal form ${\overline{H}}_n^{*}\left(J_2,{\psi }_2;J_1\right)$ is again a 1-DOF system with a parameter J₁, and ψ₂ ≈ θ₂ representing the ν₅ resonance angle. For fixed J₁, the dynamics of ${\overline{H}}_n^{*}\left(J_2,{\psi }_2;J_1\right)$ is known completely by the level curves of ${\overline{H}}_n^{*}\left(J_2,{\psi }_2;J_1\right)$, which coincide with the phase flows on the (J₂, ψ₂) plane. However, the phase plane structure now depends on the normalization order n. In this paper, we carry out the normalization algorithm to order four. Moreover, we define a simplest ideal resonance model, denoted as ${\overline{H}}_0^{*}$, by $${\overline{H}}_0^{*}\left(J_2,{\psi }_2;J_1\right)=K_0(\mathit{J})+c_{0,1}^1\cos {\psi }_2.$$(20)

Figure 3 presents the phase portraits of ${\overline{H}}_n^{*}$ with n = 0, 1,…, 4 for a particular value of J₁. We find that: the stability of the two equilibriums located at ψ₂ = 0 changes as the normalization order increases from zero to two and remains the same for n = 3, 4; (2^†). A new center at ψ₂ = π appears in frame c for n = 2, disappears in frame d for n = 3, and reappears in frame e for n = 4.

First, we recall that the level curve of J₁ is almost vertical, and the larger (smaller) the value of J₂ ≈ I₂ = P₂, the larger (smaller) the eccentricities. In fact, the level curve of J₁ adopted in Fig. 3 is located at a₁ ≈ 3.17 AU; and the upper (lower) bound of J₂ in Fig. 3 corresponds to e₁ ≈ 0.55 (e₁ ≈ 0.27). From Fig. 2, one may see that the lower bound of J₂ is touching the yellow left boundary of the J2:1. Second, we note that the closer to the libration boundary of J2:1, the slower the Fourier expansion Eq. (10) converges (see Fig. B.4). So, the phenomenon (1^†) above actually indicates the necessity of the normalization algorithm when close to the J2:1 libration boundary: the qualitative difference between frames a and b in Fig. 3 indicates that the strength of the harmonics of the form cos (lψ₂) is so strong that they cannot be simply ignored. Moreover, the phenomenon (2^†) indicates that when too close to the libration boundary, a low-order truncation of the partial normal form (such as to order four) is not convincing.

Further, it is necessary to show that when not too close to the J2:1 libration boundary, the low-order normal form is asymptotic, that is, (1*) the structures of the normal form truncated at low orders such as order four indeed reflect properly the original dynamics of ${\overline{H}}_{\text{Schubart }}$; (2*) the accuracy of the normal form increases as the normalization order increases. This part of work is tedious and saved in Appendix E, where we computed the surfaces of section of the systems $\overline{H}(\overline{P},\overline{Q})$ and $\overline{H}(I,\theta )$ to show (1*), and verify (2*) by numerical integration of example orbits.

We move on to the global structures of the ν₅ resonance inside the J2:1. For fixed J₁, we determine the locations of the equilibrium on the phase portrait and its stability for ${\overline{H}}_n^{*}$ with n = 0, 1,…,4. For centers (if exist) located at ψ₂ = ψ_2,c with ψ_2,c = 0, or, π, we determine its libration width by finding the two intersection points of the largest libration orbit with the section ψ₂ = ψ_2c. So, both the locations of centers and the libra-tion boundaries are presented by one-dimensional curves in the (J₁, J₂) plane. Moreover, instead of the (J₁, J₂) plane, it is more intuitive to present the resonance structures on the ${\overline{a}}_1-{\overline{e}}_1$ plane as in Fig. 2. To realize this, for a point (J₁, J₂) at the section ψ₂ = ψ₂,_c, we conduct the transformation⁵ $$\left({\overline{P}}_1,{\overline{P}}_2,{\overline{Q}}_1,{\overline{Q}}_2\right)={ℒ}_2°{ℒ}_{3,n}\left(J_1,J_2,{\psi }_1=0,{\psi }_{2,c}\right),$$(21)

and solve $\left({\overline{a}}_1,{\overline{e}}_1\right)$ from $\left({\overline{P}}_1,{\overline{P}}_2\right)$. The results are presented in Fig. 4.

In each frame of Fig. 4, equilibriums of the ideal resonance model ${\overline{H}}_0^{*}$ located at ψ₂ = 0 (ψ₂ = π) as well as the resonance width are presented as solid red (blue) curve as a reference for comparison. The resonance center of the normal form ${\overline{H}}_n^{*}$ with n = 1,…, 4 located at ψ₂ = 0 (ψ₂ = π) is presented as dashed red (blue) curves in frames a-d. In frame a, there exist qualitative differences for centers at ψ₂ = 0, which correspond to the phenomenon (1^†) mentioned earlier. This again indicates the non-negligible strengths of the harmonics cos (lψ₂) with l ≥ 2 in the Fourier expansion and the necessity of the normalization procedure. However, it is interesting to find that the structures of the ideal resonance model agrees well with those of the normal form ${\overline{H}}_{n\ge 2}^{*}$ as shown in frames b-d, except the structures for centers at ψ₂ = π close to the left boundary of the J2:1 resonance (the black dashed curve). As mentioned earlier by the phenomenon (2^†), the structures close to the libration boundary are suspicious because of the slow convergence (or simply divergence) of the normal form when the Fourier expansion converges slowly. In fact, results of the surfaces of section in Appendix E shows that the suspicious structures do not exist in the original systems $\overline{H}(\overline{P},\overline{Q})$ and $\overline{H}(I,\theta )$ and are thus fictitious. The boundaries of libration domain of ν₅ resonance are presented in Fig. 4 as pink (light blue) curves for centers at ψ₂ = 0 (ψ₂ = π)·We note that libration boundaries for the suspicious centers identified earlier are not presented here.

As shown in Fig. 4, the phase location of the v₅ stable centers exhibits two shifts as ${\overline{e}}_1$ increases: the first one from ψ₂ = 0 to ψ₂ = π occurs at ${\overline{e}}_1\approx 0.4$, and the second one from ψ₂ = π to ψ₂ = 0 occurs at ${\overline{e}}_1\approx 0.6$. We note that only the second one is reported in Morbidelli & Moons (1993), where in their Fig. 3, though not stated clearly, the locations of unperturbed equilibrium of v₅ resonance (namely results of K₀(J) in Eq. (20)) are presented by a smooth curve there instead of the perturbed equilibriums presented here.

Fig. 3 Phase portraits of the integrable approximation ${\overline{H}}_n^{*}\left(J_2,{\psi }_2;J_1\right)$ with n = 0, 1,…, 4 for the ν5 resonance with J1 = –1.69260419 × 10–2.

Fig. 4 Locations of v¡ resonance centers and the corresponding boundaries of the libration domain of the partial normal form ${\overline{H}}_n^{*}$ with n = 0,1,…, 4. The resonance center at ψ2 = 0 (ψ2 = π) of the ideal resonance model ${\overline{H}}_0^{*}$ is presented as a solid red (blue) curve in each frame. The boundaries of the libration domains that embrace the curves of center at Ψ2 = 0 (Ψ2 = π) are presented as pink (light blue) curves. The resonance structures of the partial normal form ${\overline{H}}_n^{*}$ with n = 1,…, 4 are presented as dashed curves in frames a-d. The vertical gray curves are level curves of I1, which approximates the level curves of J1, along which the v5 libration follows. In frame c, the yellow dot marks the pseudo-proper elements of the flow of ${\overline{H}}_{\text{Schubart }}$ in frame c of Fig. E.1.

3.2 Numerical computation of pseudo-proper elements

For an orbit of the complete system ERTBP, supposing it is trapped in the v₅ resonance, to determine its location in the libration domain in Fig. 4 is equivalent to the determination of its proper elements using the perturbation scheme developed in Sect. 2. That is, for an initial condition in terms of the chart (P, Q), we calculate $(J,\psi )={ℒ}_3^{-1}°{ℒ}_2^{-1}°{ℒ}_1^{-1}(P,Q)$; then, the proper element (J_1,pro, J_2,pro) is given by the intersection point of the flow determined by (J,ψ) with the section ψ₂ = 0, or, π. Using Eq. (21), we obtain the proper element $\left({\overline{P}}_{1,\text{pro}},{\overline{P}}_{2,\text{pro}}\right)$, from which (a_1,pro, e_1,pro) is solved, and can be projected onto Fig. 4. To conclude, the proper element is defined as the intersection of the flow of the system ${\overline{H}}_{\text{Normal }}^{*}$ with the section $${\overline{H}}_{\text{Normal }}^{*}:{\psi }_1=0\wedge {\psi }_2=0,\text{ or },\pi .$$(22)

An alternative approach is to define the so-called pseudo-proper elements (e.g., Roig et al. 2002) by numerically integrating the flow and recording intersection points using sections equivalent to (or close to) Eq. (22). To be more specific, the pseudo-proper elements can be generated using the following criteria suitable for different systems $$\{\begin{array}{cc}{H}_{\text{pERTBP }}:& Q_1\approx 0\wedge Q_2\approx 0,\text{ or },\pi \\ {\overline{H}}_{\text{Schubart }}:& {\overline{Q}}_1\approx 0\wedge {\overline{Q}}_2\approx 0,\text{ or },\pi \\ {\overline{H}}_{\text{Fourier }}:& {\theta }_1\approx 0\wedge {\theta }_2\approx 0,\text{ or },\pi \end{array}$$(23)

As a first example, the yellow dots⁶ in frame c of Fig. 4 present the pseudo-proper elements of the orbit of the system ${\overline{H}}_{\text{Schubart }}$ presented in frame c of Fig. E.1. The initial condition corresponds to a v₅ resonance center provided by the normal form of order three. As expected, the projection of the pseudo proper elements dwells exactly on the dashed blue curves instead of the solid blue curve. As noted in footnote 5, the Lie mapping ℒ₃, although a near-identity transformation, is non-negligible for an agreement between the yellow dots (pseudo-proper element) and the dashed blue curve (curve of resonance center of normal form) observed in frame c.

Next, we used the SWIFT package (Levison & Duncan 1994; Brož et al. 2011) to carry out a more systematic computation of the pseudo-proper elements. In this way, we completed the test of the first chain of the perturbation scheme in Fig. 1 by directly comparing the system ${\overline{H}}_{\text{Normal }}^{*}$. First, we generated a grid $S_3=\left\{\left({\overline{a}}_{1,(i,j)},{\overline{e}}_{1,(i,j)}\right)\mid 1\le i\le 10,1\le j\le 20\right\}$ that covers the libration domain of the J2:1 as shown by black crosses in Fig. 5. Using the relations in Eq. (2), we derive the corresponding $\left({\overline{P}}_{1,(i,j)},{\overline{P}}_{2,(i,j)}\right)$ from $\left({\overline{a}}_{1,(i,j)},{\overline{e}}_{1,(i,j)}\right)$. Then, fixing that ${\overline{Q}}_1=0$ and ${\overline{Q}}_2=0$, or, π, we obtain the corresponding (P, Q)_(i,j) from ${(\overline{P},\overline{Q})}_{(i,j)}$ using ℒ₁ under the condition that Q₃ = Q₄ = 0. In this way, we obtained two sets of initial conditions that both cover the libration domain of the J2:1. One set satisfies ɡ₁ − ɡ₂ ≈ 0 and the other set ɡ₁ − ɡ₂ ≈ π. Then, the two sets of initial conditions are numerically integrated using the SWIFT package for 1 Myr to record its pseudo-proper elements using the first criterion for H_pERTBP in Eq. (23). The results are presented as dots in Fig. 5. As shown in Fig. 5, for orbits whose initial conditions are outside the ν₅ libration domain, the gray dots generally coincide with the black crosses as expected. For orbits initially dwelling in the ν5 libration domain (colored dots), there are two clusters of pseudo-proper elements that are distributed on the two sides of the curve of ν₅ resonance center along the level curves of I₁ ; and one of the two clusters coincide with the black crosses. The initial conditions that have no results of pseudo-proper elements correspond to very unstable orbits that escape the J2:1 within a short-time interval.

Fig. 5 Pseudo-proper elements (dots) of orbits inside the J2:1 commen-surability obtained by the SWIFT package. The black crosses mark their initial conditions with 𝑔1 – 𝑔2 ≈ 0 (𝑔1 – 𝑔2 = π) in the left (right) frame. The section Q2 ≈ 0 (Q2 ≈ π) in Eq. (23) is used to record the pseudo-proper elements of each orbit in frame a (frame b). The colored (gray) dots mark the pseudo-proper elements of orbits trapped (not trapped) in the v5 libration domain.

Fig. 6 Drifts of pseudo-proper elements of orbits inside the J2:1 commensurability due to the diurnal Yarkovsky effect. The initial conditions are the same as those in Fig. 5 while the Yarkovsky effect is included in the force model (see the text for physical parameters related to the Yarkovsky effect). As in Fig. 5, the section Q2 ≈ 0 (Q2 ≈ π) in Eq. (23) is used to record the pseudo-proper elements of each orbit in frame a (frame b). The integration time is 2 Myr. The ν5 resonance structures are taken from Fig. 4. The gray curves are the level curves of I1 that the ν5 libration approximately follows. Compared with Fig. 5, the pseudo-proper elements show obvious systematic drifts. The light yellow curves in frame a mark the pseudo-proper elements obtained using partial normal form H2* perturbed by ${{\overline{F}}_Y|}_J$, which agrees well with the reference results by the SWIFT package.

4 Drifts in the action space of the J2:1 resonance by the Yarkovsky effect

4.1 Classification of systematic drifts by numerical simulations

We started by running numerical simulations using the two sets of initial conditions in Fig. 5 and the augmented version of the SWIFT package by Brož et al. (2011), which supports the option of the Yarkovsky effect in the force model. The physical parameters for the Yarkovsky effect are chosen as R = 0.5 m, ω = ω(R = 0.5 m), ρ = ρ_reg, K = 10⁻⁶ W/ (m K), C = C_rock (see Eq. (D.5)), and θ₀ = 0, that is, all asteroids are assumed to be prograde rotators. Then, we record the pseudo-proper elements of each orbit during 2 Myr integration (using the section for H_pERTBP in Eq. (23)) and present the results in Fig. 6 in the same pattern as in Fig. 5.

Compared with Fig. 5, it is clear that the pseudo-proper elements show systematic drifts because of the Yarkovsky effect. The drift patterns can be generally classified into four types:

• Type Ia (Light gray dots): the orbits are not trapped in the ν5 resonance and the drifts are almost parallel to the level curves of I₁ ;

• Type Ib (Dark gray dots): the orbits are not trapped in the ν₅ resonance and the drifts cross the level curves of I₁ ;

• Type IIa (colored dots except red/blue): the orbits are trapped in the ν₅ resonance and the drifts are almost parallel to the ν5 resonance curves;

• Type IIb (red and blue dots): the orbits are temporarily trapped in the ν₅ resonance during the integration and the drifts are almost parallel to the level curves of I₁ but “jumps” over the ν₅ libration domain.

To better understand the above classifications, Fig. 7 presents for each type an example trajectory on the plane (ɡ₁ − ɡ₂ , e₁). Frame a presents two orbits of Type Ia marked as gray dots with black edges in the left frame of Fig. 6: the upper (lower) panel corresponds to the orbit with ${\overline{a}}_1=3.19\text{\hspace{0.17em}AU\hspace{0.17em}}\left({\overline{a}}_1=3.26\text{\hspace{0.17em}AU\hspace{0.17em}}\right)$. The ν5 resonance angle ɡ₁ − ɡ₂ is in circulation with d/dt (ɡ₁ − ɡ₂) < 0 for both orbits while the peak of eccentricity is located at ɡ₁ − ɡ₂ = π and ɡ₁ − ɡ₂ = 0 for the upper panel and the lower panel, respectively (The reason is easily explained using the ν₅ resonance portrait for the corresponding I₁ value of each orbit). For this reason, the orbits with larger J2:1 resonance libration amplitudes appear to drift faster because the section to record pseudo-proper elements is chosen as ɡ₁ − ɡ₂ = 0 in the left frame of Fig. 6. The same phenomenon exists in the right frame of Fig. 6 as shown by the two orbits marked by gray dots with edges. However, the situation is converse now: the orbits with smaller libration amplitude appear to drift faster because the section to record pseudo-proper elements is chosen as ɡ₁ − ɡ₂ = π. We may conclude that for Type Ia drifts, the pseudo-proper semi-major axis does not change much, and the drift rate of the eccentricity has a weak dependence on the J2:1 resonance libration amplitude.

Frame b presents two orbits of Type IIa where the top (bottom) panel corresponds to the yellow (brown) dots in the libration domain of Ψ₂ = 0 (Ψ₂ = π) with ${\overline{a}}_1=3.17\text{\hspace{0.17em}AU\hspace{0.17em}}\left({\overline{a}}_1=3.18\text{\hspace{0.17em}AU\hspace{0.17em}}\right)$ in the left (right) frame of Fig. 6. It is clear that ɡ₁ − ɡ₂ is in libration for both orbits. The difference is that in the upper panel the libration cycle is going up while in the lower panel the libration cycle is going down as marked by the arrows. In fact, the common characteristic for drifts of Type IIa in Fig. 5 is that they all move in the leftward direction in the pseudo-proper a₁ − e₁ plane with the footprints almost parallel to the v₅ resonance curve. So, for an orbit trapped in the v₅ libration domain with centers located at Ψ2 ≈ g₁ − g₂ = 0 (namely the regions encircled by the pink curves) at about 3.17 AU, its eccentricity would increase as it drifts parallel to the resonance curve in the leftward direction. On the other hand, for an orbit trapped in the libration domain encircled by pink curves at about 3.25 AU, its eccentricity would decrease. For orbits trapped in the v₅ libration domain encircled by the light blue curves, its eccentricity would always decrease.

Frame c presents an orbit of type Ib, which corresponds to the dark gray dots with black edges located at the top of the left frame of Fig. 6. From frame a of Fig. 6, the orbit approaches the stable libration boundary of J2:1 as the pseudo-proper semimajor axis decreases due to the Yarkovsky force while the eccentricity does not show obvious growth. Frame c shows that g₁ − g₂ is in circulation with d/dt (g₁ − g₂) > 0 under the Yarkovsky force. As the orbit hits the boundary, the close encounter with Jupiter changes its semi-major axis and eccentricity abruptly. Then, the eccentricity undergoes large-amplitude oscillation due to secular resonances.

Frame d presents two orbits of Type IIb where the top (bottom) panel corresponds to red (blue) dots at 3.25 AU (3.18 AU) in the left frame of Fig. 5. Both orbits start at the section g₁ − g₂ = 0, and its evolution exhibits three phases: The first phase is marked by the arrow 1, during which g₁ − g₂ is in circulation with d/dt (g₁ − g₂) < 0 and the eccentricity increases slowly, namely the Type Ia drift. The second phase is marked by the arrow 2, during which the orbit is temporarily trapped in v₅ resonance, and the eccentricity increases rapidly during one incomplete cycle of the v₅ libration. The last phase, marked by arrow 3, is again a Type Ia drift in which the g₁ − g₂ is in circulation with d/dt (g₁ − g₂) > 0 and the eccentricity increases slowly. Intuitively speaking, the orbit of Type IIb jumps over the libration domain under the Yarkovsky effect and gains a rapid increase of eccentricity as a result.

Fig. 7 Trajectories on the (e 1,ɡ1 − ɡ2) plane for example orbits of different types of drift identified in Fig. 6. Frame a presents two orbits of Type Ia marked as gray dots with black edges in the left frame of Fig. 6; Frame b presents two orbits of Type IIa dwelling in the libration domain of Ψ2 = 0 (Ψ2 = π) at ${\overline{a}}_1=3.17\text{\hspace{0.17em}AU\hspace{0.17em}}\left({\overline{a}}_1=3.18\text{\hspace{0.17em}AU\hspace{0.17em}}\right)$ in the left (right) frame of Fig. 6. Frame c presents an orbit of Type Ib, marked as dark gray dots with black edges in the left frame of Fig. 6. Frame d presents two orbits of Type IIb where the top (bottom) panel corresponds to red (blue) dots at ${\overline{a}}_1=3.17\text{\hspace{0.17em}AU\hspace{0.17em}}\left({\overline{a}}_1=3.18\text{\hspace{0.17em}AU\hspace{0.17em}}\right)$ in the left frame of Fig. 6. In frames b and d, the arrows mark the directions of evolution.

4.2 Analysis on the drift rates of pseudo-proper elements

We compute in Sect. 2.5 the averaged values of the differentiate rates of J₁, J₂ due to the Yarkovsky force, namely ${\overline{F}}_Y|\mathit{J}$ . Our goal is to link ${\overline{F}}_{{Y|}_{\mathit{J}}}$ with the drift rates of the pseudo-proper semi-major axis and eccentricity, namely $\text{d}{a}_{1,\text{pse}}\text{/d}t\text{\hspace{0.17em}and\hspace{0.17em}d}{e}_{1,\text{pse}}\text{/d}t$ and explain the different types of drifts observed in the previous section.

Before moving on, it is necessary to first show that the several types of drifts identified above indeed exist in the nor- malized systems ${\overline{H}}_n^{*}(J,\psi )$ under the perturbation of the averaged Yarkovsky force ${{\overline{F}}_Y|}_{\mathit{J}}$. So, we integrate the initial conditions of several orbits in Fig. 6 using the canonical equations of H₂^* perturbed by ${\overline{F}}_{{Y|}_{\mathit{J}}}$ and record the pseudo-proper elements, which are presented in the left frame of Fig. 7 using small light yellow dots (which appear as thin yellow curves). It is shown that the yellow dots agree well with the corresponding colored dots. This suggests that the dynamics that interests us is preserved in the normalized system with the averaged Yarkovsky force.

With Eq. (21), we defined two tabulated functions a_1,pse (J_1,pse, J_2,pse), and e_1,pse J_1,pse, J_2,pse) on a grid in the action space J₁ − J₂ . Then, the differentiation of the two functions bridges the drift rates of (a_1,pse, e_1,pse) and that of J_1,pse, J_2,pse). Because J₁,_pse is a parameter of the partial normal form ${\overline{H}}_n^{*}$, we have J_1,pse = J₁ and thus $\text{d}{J}_{1,\text{pse}}/\text{d}t={\overline{F}}_Y{|}_{J_1}$.

We recall by definition that J_2,pse satisfies the relation H(J_2,pse,Ψ₂ = 0; J₁) = h. So we write J_2,pse = J_2,pse(J₁, h), for which the differentiation gives (similar to Eq. (16)) $$\begin{array}{l}\frac{\text{d}{J}_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\text{d}t}={\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial J_1}{\overline{F}}_Y|}_{J_1}\\ +\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial h}\left(\frac{\partial H_n^{*}}{\partial J_1}{\overline{F}}_{{Y|}_{J_1}}+\frac{\partial H_n^{*}}{\partial J_2}{\overline{F}}_{{Y|}_{J_2}}+\frac{\partial H_n^{*}}{\partial {\psi }_2}{\overline{F}}_{Y{\psi }_{{\psi }_2}}\right).\end{array}$$(24)

For applications, we should average Eq. (24) over one cycle (libration or circulation) of the v₅ resonance, namely $\langle \text{d}{J}_{2,\text{pse}}/\text{d}t\rangle =1/T_{v_5}{\displaystyle \oint \left(\text{d}{J}_{2,\text{pse}}/\text{d}t\right)}\text{d}t$, where T_v5 is the period of the v₅ resonance. We note that the two cases of the integration trajectory, namely circulation and libration, correspond to the Type I and the Type II drifts, respectively.

Since ${\overline{F}}_Y$ is very small, the averaging can be carried along a flow of the normal form without considering the effects of ${\overline{F}}_Y$ . That is, J1 and h are assumed as constants in the right-hand side of Eq. (24). So, the variations of ${\overline{F}}_Y$ during one cycle of v₅ resonance is subject only to the oscillation of J₂ . The contour maps of $\left|{\overline{\mathit{F}}}_Y{|}_J\right|$ (not shown here for simplicity) shows that (1) the variation of ${\overline{F}}_Y{|}_{J_1}$ is small for fixed a1 (which can approximate the level curve of J₁ along which the v₅ libration follows); (2) the variation of ${{\overline{F}}_Y|}_{J_2}$ is small in the whole J2:1 libration domain. So, it is reasonable to assume that ${{\overline{F}}_Y|}_J$ is constant during one cycle of v₅ resonance so that $$\begin{array}{l}\langle \frac{\text{d}{J}_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\text{d}t}\rangle =\left(\frac{\partial J_{2,\text{pse}}}{\partial J_1}+\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial h}\langle \frac{\partial H_n^{*}}{\partial J_1}\rangle \right){\overline{F}}_{{Y|}_{J_1}}\\ +\frac{\partial J_{2,\text{pse}}}{\partial h}\left(\langle \frac{\partial H_n^{*}}{\partial J_2}\rangle {\overline{F}}_{{Y|}_{J_2}}+\langle \frac{\partial H_n^{*}}{\partial {\psi }_2}\rangle {\overline{F}}_{{Y|}_{{\psi }_2}}\right).\end{array}$$(25)

4.2.1 Type Ia drift

Consider again the relation H(J_2,pse,Ψ₂=0;J₁)=h. Taking J_2,pse as a function of h and J₁ , the differentiation of the formula gives $$\{\begin{array}{l}\frac{\partial H}{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial J_1}+\frac{\partial H}{\partial J_1}=0\hfill \\ \frac{\partial H}{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial h}=1\hfill \end{array}\Rightarrow \frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial J_1}+\frac{\partial J_{2,\text{\hspace{0.17em}pse\hspace{0.17em}}}}{\partial h}\frac{\partial H}{\partial J_1}=0.$$(26)