An analytically treatable model of long-term dynamics in a mean motion resonance with coexisting resonant modes

Abstract

We examine a 2-DOF Hamiltonian system modeling some previously unexplored dynamical effects in first-order mean motion resonances in the spatial circular restricted three-body problem “star-planet-asteroid.” In distinction from the well-known integrable model of Sessin and Ferraz-Mello (Celest Mech Dyn Astron 32:307–332, 1984) and Wisdom (Celest Mech Dyn Astron 38:175–180, 1986), in our analysis, we kept more terms in the expansion for the disturbing function. This allowed us to study the phenomena caused by coexisting resonant modes, that were lost in the integrable model. In addition, we identified the areas of dynamical chaos in the phase space of the considered system. Despite the non-integrability of our model, we were able to compute analytically many quantities characterizing its dynamical behavior. Finally, we illustrate some of our results through the examples of the known resonant Kuiper belt objects.

Appendices

Appendix A: Constructing a model system, that reveals the origin of chaos in first-order MMR

We shall confine ourselves to the case of an exterior resonance $p:(p+1)$ within the restricted three-body problem. The interior resonance $(p+1):p$ can be reduced to the same model by using similar approach. The distance between the two major bodies (i.e. a star and a planet) and the sum of their masses are taken here as units of length and mass. The unit of time is chosen such that the orbital period of the major bodies’ rotation about their barycenter is equal to $2\pi$. The mass of the planet $\mu$ is considered to be a small parameter of the problem.

The equations of motion for the minor body (asteroid) in the canonical form are

$$\begin{aligned} \frac{\mathrm{d}(L,G,H)}{\mathrm{d}t} = - \frac{\partial \mathcal {K}}{\partial (l,g,h)},\quad \frac{\mathrm{d}(l,g,h)}{\mathrm{d}t} = \frac{\partial \mathcal {K}}{\partial (L,G,H)}, \end{aligned}$$

(34)

where L, G, H, l, g, h are the Delaunay variables (Murray and Dermott 2000). They can be expressed through the Keplerian elements $a,e,i,\varOmega , \omega$ as

$$\begin{aligned} L = \sqrt{(1 - \mu )a},\quad G = L\sqrt{1 - {e^2}},\quad H = G\cos i,\quad g = \omega ,\quad h = \varOmega . \end{aligned}$$

The last variable l is the mean anomaly of the asteroid.

The Hamiltonian $\mathcal {K}$ in (34) is

$$\begin{aligned} \mathcal {K} = - \frac{{{{(1 - \mu )}^2}}}{{2{L^2}}} - \mu R(L,G,H,l,g,h - \lambda ^\prime ). \end{aligned}$$

(35)

Here, R is the disturbing function in the restricted circular three-body problem. The mean longitude of the planet $\lambda ^\prime$ appearing in (35) is linear in time: $\lambda ^\prime = t + {\lambda ^\prime _0}$. Therefore, it is convenient to employ the variable ${\tilde{h}} = h - \lambda ^\prime$ instead of h, as it enables us to write down the equations of motion in an autonomous form as canonical equations with the Hamiltonian

$$\begin{aligned} \tilde{ \mathcal {K}} = \mathcal {K}(L,G,H,l,g,{\tilde{h}}) - H. \end{aligned}$$

We introduce a resonant angle ${\bar{\varphi }}$ by using the canonical transformation $(L,G,H,l,g,{\tilde{h}}) \rightarrow ({P_\varphi },{P_g},{P_h},{\bar{\varphi }} ,{\bar{g}},{\bar{h}})$ defined by the generating function

$$\begin{aligned} S = (p + 1){P_\varphi }l + \left[ {{P_h} + p\left( {{P_\varphi } - P_\varphi ^*} \right) } \right] {\tilde{h}} + \left[ {{P_g} + (p + 1)\left( {{P_\varphi } - P_\varphi ^*} \right) } \right] g, \end{aligned}$$

where $P_\varphi ^* = L_*/(p+1)$, while $L_*=\root 3 \of {(p+1)/p}$ is the value of L corresponding to an exact $p:(p+1)$ MMR in the unperturbed problem ($\mu =0$). The new variables are related to the old ones via

$$\begin{aligned}&\displaystyle L = \frac{{\partial S}}{{\partial l}} = (p + 1){P_\varphi },\quad G = \frac{{\partial S}}{{\partial g}} = {P_g} + (p + 1)\left( {{P_\varphi } - P_\varphi ^*} \right) ,\\&\displaystyle H = \frac{{\partial S}}{{\partial {\tilde{h}}}} = {P_h} + p\left( {{P_\varphi } - P_\varphi ^*} \right) ,\\&\displaystyle {\bar{\varphi }} = \frac{{\partial S}}{{\partial {P_\varphi }}} = (p + 1)l + p{\tilde{h}} + (p + 1)g,\quad {\bar{g}} = \frac{{\partial S}}{{\partial {P_g}}} = g,\quad {\bar{h}} = \frac{{\partial S}}{{\partial {P_h}}} = {\tilde{h}}. \end{aligned}$$

The resonant angle can also be expressed as

$$\begin{aligned} {\bar{\varphi }} = (p + 1)\lambda - p\lambda ^\prime - \varOmega . \end{aligned}$$

Although this definition does not satisfy D’Alembert rule on the sum of coefficients before longitudes of ascending nodes, using introduced variable instead of the traditional resonant angle results in a simpler form of the final system, while the averaging methods we use to study the system are applicable just the same. The new Hamiltonian is:

$$\begin{aligned} \tilde{{\mathcal {K}}}=-\frac{(1-\mu )^2}{2(p+1)^2 P_\varphi ^2}-\left[ P_h+p(P_\varphi -P_\varphi ^*)\right] - \mu R. \end{aligned}$$

The resonant case we are interested in corresponds to a region $\mathcal {R}$ of phase space, defined by the condition

$$\begin{aligned} \left| {pn^\prime - (p + 1)n} \right| \lesssim {\mu ^{1/2}}. \end{aligned}$$

Here, $n^\prime =1$ and n are the mean motions of the planet and asteroid, respectively. It is also true in $\mathcal {R}$ that

$$\begin{aligned} \left| {{P_\varphi } - P_\varphi ^*} \right| \lesssim {\mu ^{1/2}}. \end{aligned}$$

(36)

In the resonant case, i.e., when the previous inequality holds, the variables can be divided into fast, semi-fast, and slow. The fast and semi-fast variables in $\mathcal {R}$ are ${\bar{h}}$ and ${\bar{\varphi }}$ respectively:

$$\begin{aligned} \frac{{\mathrm{d}{\bar{h}}}}{{\mathrm{d}t}} \sim 1,\quad \frac{{\mathrm{d}{\bar{\varphi }} }}{{\mathrm{d}t}} \sim {\mu ^{1/2}}. \end{aligned}$$

The slow variables (whose rate is of order $\mu$), are $P_\varphi$, $P_g$, $P_h$, and ${\bar{g}}$.

To study secular effects, averaging over the fast variable ${\bar{h}}$ is performed, which results in equations of motion assuming a canonical form with the Hamiltonian

$$\begin{aligned}&{\bar{ {\mathcal {K}}}}({P_\varphi },{P_g},{P_h},{\bar{\varphi }} ,{\bar{g}}) \nonumber \\&\quad =\frac{1}{{2\pi (p + 1)}}\int \limits _0^{2\pi (p + 1)} {\tilde{{\mathcal {K}}}(L({P_\varphi }),G({P_g},{P_\varphi }),H({P_h},{P_\varphi }),l({\bar{\varphi }} ,{\bar{g}},{\bar{h}}),{\bar{g}},\bar{h})} \,\mathrm{d}{\bar{h}}, \end{aligned}$$

where

$$\begin{aligned} l({\bar{\varphi }} ,{\bar{g}},{\bar{h}}) = \frac{1}{{p + 1}}\left[ {{\bar{\varphi }} - (p + 1){\bar{g}} - p{\bar{h}}} \right] . \end{aligned}$$

After such averaging, the fast variable $\bar{h}$ vanishes, and the term “fast” becomes redundant. Therefore in the rest of the paper we adopt the name fast to denote the variables which vary with the rate $\mu ^{1/2}$, instead of referring to them as semi-fast (as we used to do when the distinction among three different time scales was needed).

Moreover, since ${\bar{h}}$ is no longer present in the expression for the Hamiltonian $\bar{ {\mathcal {K}}}$, the conjugate momentum $P_h$ is a constant in the considered approximation and can be treated as a parameter of the problem. Thus, $\bar{ {\mathcal {K}}}$ is the Hamiltonian of a system with two degrees of freedom. Further on, instead of $P_h$ we shall use the parameter

$$\begin{aligned} \sigma = \sqrt{1 - \frac{P_h^2}{L_*^2}}. \end{aligned}$$

The inequality $e\le \sigma$ defines the region ${\mathcal {S}}$ in the phase space, to which the motion of the system is bound by the Lidov–Kozai integral (Sidorenko et al. 2014).

The next standard step in the analysis of the dynamics in the resonant region ${\mathcal {R}}$ is the scaling transformation (Arnold et al. 2006):

$$\begin{aligned} {\bar{\tau }} = {\bar{\varepsilon }} t,\quad {\bar{\varPhi }} = (P_\varphi ^* - {P_\varphi })/{\bar{\varepsilon }}, \end{aligned}$$

(37)

where ${\bar{\varepsilon }} = {\mu ^{1/2}}$ is a new small parameter. Through using variables (37), the equations of motion can be rewritten, without loss of accuracy, in the form of a slow-fast system:

$$\begin{aligned} \frac{\mathrm{d}{\bar{\varphi }}}{\mathrm{d}{\bar{\tau }}}= & {} \chi {\bar{\varPhi }},\quad \frac{\mathrm{d}{\bar{\varPhi }}}{\mathrm{d}{\bar{\tau }}} = - \frac{{\partial \bar{W}}}{{\partial {\bar{\varphi }} }},\nonumber \\ \frac{\mathrm{d} P_g}{\mathrm{d}{\bar{\tau }}}= & {} {\bar{\varepsilon }} \frac{{\partial \bar{W}}}{{\partial {\bar{g}}}},\quad \frac{\mathrm{d}{\bar{g}}}{\mathrm{d}{\bar{\tau }}} = - {\bar{\varepsilon }} \frac{{\partial \bar{W}}}{{\partial {P_g}}}, \end{aligned}$$

(38)

where

$$\begin{aligned} \begin{aligned}&\chi = 3{p^{4/3}}{\left( p + 1\right) ^{2/3}}, \\&\qquad \bar{W}({\bar{\varphi }} ,{\bar{g}},{P_g};\sigma ) = \frac{1}{{2\pi (p + 1)}} \int \limits _0^{2\pi (p + 1)} {R\left( {L_*},{P_g},{L_*}\sqrt{1 - {\sigma ^2}} ,l({\bar{\varphi }} ,{\bar{g}},{\bar{h}}),{\bar{g}},{\bar{h}} \right) } \,\mathrm{d}{\bar{h}}. \end{aligned} \end{aligned}$$

For $\sigma \ll 1$, an approximate expression for $\bar{W}(\bar{\varphi },{\bar{g}},{P_g};\sigma )$ can be written down as a series expansion (Murray and Dermott 2000):

$$\begin{aligned}&\bar{W}({\bar{\varphi }} ,{\bar{g}},{P_g};\sigma ) \approx {W_0} + {W_1}e\cos ({\bar{\varphi }} - {\bar{g}}) \nonumber \\&\quad + {e^2}\left[ {{W_{10}} + {W_{11}}\cos 2({\bar{\varphi }} - {\bar{g}})} \right] + {i^2}\left[ {W_{20} + {W_{21}}\cos 2{\bar{\varphi }} } \right] , \end{aligned}$$

(39)

where

$$\begin{aligned} e^2\approx 1-\frac{P_g^2}{L_*^2},\quad i^2\approx 2\left( 1-\frac{P_h}{L_* \sqrt{1-e^2}}\right) \approx \sigma ^2-e^2. \end{aligned}$$

(40)

Expressions (40) are obtained by using the resonance condition (36) and by recalling that the values e and i are limited by $\sigma \ll 1$, leading to $e, i\ll 1$.

The coefficients $W_0$, $W_1$, $W_{10}$, $W_{11}$, $W_{20}$, $W_{21}$ in (39) are calculated as follows:

$$\begin{aligned} \begin{array}{l} W_0 = \frac{\alpha }{2}b_{1/2}^{(0)},\quad W_1 = \frac{\alpha }{2}\left[ {(2p + 1)b_{1/2}^{(p)} + \alpha \frac{{\mathrm{d}b_{1/2}^{(p)}}}{{\mathrm{d}\alpha }}}\right] - \frac{\delta _{1p}}{2\alpha },\\ W_{10} = \frac{\alpha ^2}{8}\left( {2\frac{{\mathrm{d}b_{1/2}^{(0)}}}{{\mathrm{d}\alpha }} + \alpha \frac{{{\mathrm{d}^2}b_{1/2}^{(0)}}}{{\mathrm{d}{\alpha ^2}}}} \right) ,\\ W_{11} = \frac{\alpha }{8}\left( {\left[ {2 - 14(p + 1) + 16{{(p + 1)}^2}} \right] b_{1/2}^{(2p)} + \alpha \left[ {8(p + 1) - 2} \right] \frac{{\mathrm{d}b_{1/2}^{(2p)}}}{{\mathrm{d}\alpha }} + {\alpha ^2}\frac{{{\mathrm{d}^2}b_{1/2}^{(2p)}}}{\mathrm{d}\alpha ^2}} \right) ,\\ W_{20}= - W_{10},\quad {W_{21}} = \frac{{{\alpha ^2}}}{8}b_{3/2}^{(2p + 1)}, \end{array} \end{aligned}$$

(41)

$\delta _{mn}$ is the Kronecker delta, and $b_{1/2}^{(n)}(\alpha )$, $b_{3/2}^{(n)}(\alpha )$ are the Laplace coefficients, while $\alpha = {(p/(p + 1))^{2/3}}$. The numerical values of the coefficients (41) for several resonances are gathered in Table 2.

Table 2 Numerical values of coefficients in the series expansion of the averaged disturbing function for semimajor axes ratio defined by the resonance degree p

Out of the whole region ${\mathcal {S}}$, we are interested in the part with small eccentricities. Thus, we further assume $e/\sigma \lesssim \sigma$, which leads to $e\lesssim \sigma ^2$. Taking into account (40), we can write down an expression for the averaged disturbing function up to the terms of order $\sigma ^2$:

$$\begin{aligned} \bar{W}({\bar{\varphi }} ,{\bar{g}},{P_g};\sigma ) \approx {W_0} + {W_1}e\cos ({\bar{\varphi }} - {\bar{g}}) + {\sigma ^2}\left[ {{W_{20}} + {W_{21}}\cos 2{\bar{\varphi }} } \right] . \end{aligned}$$

By introducing the variables

$$\begin{aligned} {\bar{x}} = \sqrt{2({L_*} - {P_g})} \cos {\bar{g}},\quad {\bar{y}} = - \sqrt{2({L_*} - {P_g})} \sin {\bar{g}}, \end{aligned}$$

the equations of motion (38) can be reduced to the Hamiltonian form:

$$\begin{aligned} \frac{{\mathrm{d}{\bar{\varphi }} }}{{\mathrm{d}{\bar{\tau }} }}= & {} \frac{\partial \bar{\varXi }}{{\partial {\bar{\varPhi }} }},\quad \frac{{\mathrm{d}{\bar{\varPhi }} }}{{\mathrm{d}{\bar{\tau }} }} = - \frac{\partial \bar{\varXi }}{{\partial {\bar{\varphi }} }},\\ \frac{{\mathrm{d}{\bar{x}}}}{{\mathrm{d}{\bar{\tau }} }}= & {} {\bar{\varepsilon }} \frac{\partial \bar{\varXi }}{{\partial {\bar{y}}}},\quad \frac{{\mathrm{d}{\bar{y}}}}{{\mathrm{d}{\bar{\tau }} }} = - {\bar{\varepsilon }} \frac{\partial \bar{\varXi }}{{\partial {\bar{x}}}} \end{aligned}$$

with the Hamiltonian

$$\begin{aligned} \bar{\varXi } = \frac{{\chi {\bar{\varPhi } ^2}}}{2} + \frac{{{W_1}}}{{\sqrt{{L_*}} }}\cos {\bar{\varphi }} \cdot {\bar{x}} - \frac{{{W_1}}}{{\sqrt{{L_*}} }}\sin {\bar{\varphi }} \cdot {\bar{y}} + {\sigma ^2}{W_{21}}\cos 2{\bar{\varphi }}. \end{aligned}$$

The final rescaling of the variables,

$$\begin{aligned}&\displaystyle \varphi = - \bar{\varphi },\quad \varPhi = - \sqrt{\frac{\chi }{{{\sigma ^2}{W_{21}}}}} \bar{\varPhi },\nonumber \\&\displaystyle x = \frac{{{W_1}}}{{{\sigma ^2}{W_{21}}\sqrt{{L_*}} }}\bar{x},\quad y= \frac{{{W_1}}}{{{\sigma ^2}{W_{21}}\sqrt{{L_*}} }}\bar{y},\nonumber \\&\displaystyle \tau = \sigma \sqrt{\chi {W_{21}}} \bar{\tau },\quad \varepsilon = \frac{{{W_1}^2}}{{{{\chi ^{1/2} {L_*}}}{\sigma ^3}W_{21}^{3/2}}}\bar{\varepsilon }, \end{aligned}$$

(42)

results in the model system (3)–(4).

Appendix B: Analytical expressions for the integrals emerging in the averaging over the fast subsystem’s period

The right-hand side parts of the evolution equations (16), as well as the adiabatic invariant formula (26), can be expressed in terms of the integrals

$$\begin{aligned} I_{k,r} = \int \limits _{\lambda _*}^{\lambda ^*} \frac{{\lambda ^r}\mathrm{d}\lambda }{{({\lambda ^2} + 1)^k}\sqrt{\pm (\lambda -a_1)(\lambda -a_2)(\lambda -a_3)(\lambda -a_4)}}, \end{aligned}$$

where $k = 0,1,2$; $r = 0,1$, while the integration limits $\lambda _*$, $\lambda ^*$ can be either real numbers or $\pm \infty$. The infinite limits appear when averaging is performed along librating solutions of the fast subsystem, which violate the condition (21). In this case, after the substitution $\lambda =\tan \left( \varphi /2\right)$ the integration in (22) and (26) is carried over two semi-infinite intervals. E.g., for $\varphi _*<\pi ,~ \varphi ^*>\pi ,~ \varphi ^*-\varphi _*<2\pi$, the expression for the period T becomes

$$\begin{aligned} T(x,y,\xi ) = 4\left( {\int \limits _{-\infty }^{\lambda ^*} {\frac{{\mathrm{d}\lambda }}{{\sqrt{2{R_4}(\lambda )} }}} + \int \limits _{\lambda _*}^{+\infty } {\frac{{\mathrm{d}\lambda }}{{\sqrt{2{R_4}(\lambda )} }}} } \right) . \end{aligned}$$

In such cases, it is implied that all $I_{k,r}$ on the right-hand sides of (23) and (26) are sums of two integrals as well.

The integrals $I_{k,r}$ can be reduced to the linear combinations of elliptic integrals of the first, second, and third kinds:

$$\begin{aligned} I_{0,0}= & {} c_{0,0}K(m),\nonumber \\ I_{1,0}= & {} c_{1,1}K(m) + c_{1,3}\varPi \left( {n\left| m\right. } \right) + {{{\bar{c}}}_{1,3}}\varPi \left( {{\bar{n}}\left| m\right. } \right) ,\nonumber \\ I_{1,1}= & {} g_{1,1}K(m) + g_{1,3}\varPi \left( {n\left| m\right. } \right) + {{{\bar{g}}}_{1,3}}\varPi \left( {{\bar{n}}\left| m\right. } \right) ,\nonumber \\ I_{2,0}= & {} \frac{1}{2}{I_{1,0}} + {c_{2,1}}K(m) + {c_{2,2}}E(m) + {c_{2,3}}\varPi \left( {n\left| m\right. } \right) + {{{\bar{c}}}_{2,3}}\varPi \left( {{\bar{n}}\left| m\right. } \right) ,\nonumber \\ I_{2,1}= & {} {g_{2,1}}K(m) + g_{2,2}E(m) + {g_{2,3}}\varPi \left( {n\left| m\right. } \right) + {{{\bar{g}}}_{2,3}}\varPi \left( {{\bar{n}}\left| m\right. } \right) , \end{aligned}$$

(43)

where the following notation is used:

$$\begin{aligned} K(m)= & {} \int \limits _{0}^{\pi /2} \frac{\mathrm{d}\theta }{\sqrt{1-m \sin ^2 \theta }}, \quad E(m)=\int \limits _{0}^{\pi /2} \sqrt{1-m \sin ^2 \theta }\mathrm{d}\theta , \\ \varPi \left( n\left| m\right. \right)= & {} \int \limits _{0}^{\pi /2} \frac{\mathrm{d}\theta }{\left( 1-n \sin ^2 \theta \right) \sqrt{1-m \sin ^2 \theta }}. \end{aligned}$$

Further, the formulae for coefficients $c_{k,l}$, $g_{k,l}$, parameter m and characteristic n are gathered for all necessary cases.

1.1 The case ${a_j} \in {{\mathbb {R}}^1}$ $(j = \overline{1,4} )$

We will further assume that $a_j$ are numbered in the ascending order:

$$\begin{aligned} a_1<a_2<a_3<a_4. \end{aligned}$$

Four different instances should be considered:

• A. Integration over $(a_1,a_2)$;

• B. Integration over $(a_2,a_3)$;

• C. Integration over $(a_3,a_4)$;

• D. Integration over $(-\infty ,a_1)\bigcup (a_4,+\infty )$.

1.1.1 Instance A

Here, $\lambda _*=a_1$ and $\lambda ^*=a_2$. We also shall use the following auxiliary parameters:

$$\begin{aligned} {A_0} = \frac{2}{{\sqrt{({a_4} - {a_2})({a_3} - {a_1})} }},\quad {C_0} = \frac{1}{{{a_2} - i}},\quad {\alpha ^2} = \frac{{{a_2} - {a_1}}}{{{a_3} - {a_1}}},\quad \alpha _1^2 = \frac{{({a_2} - {a_1})(i - {a_3})}}{{({a_3} - {a_1})(i - {a_2})}}. \end{aligned}$$

The value of the elliptic integrals’ parameter in (43) is:

$$\begin{aligned} m=\frac{(a_4-a_3)(a_2-a_1)}{(a_4-a_2)(a_3-a_1)}, \end{aligned}$$

while the characteristic is given by $n = \alpha _1^2$.

In order to find the coefficients in (43), we consider such linear combinations of ${I_{k,r}}$ which can be reduced to integral (253.39) from Byrd and Friedman (1954):

$$\begin{aligned} I_{1,1} + i I_{1,0}= & {} \int \limits _{a_1}^{a_2} {\frac{{\mathrm{d}\lambda }}{{(\lambda -i)\sqrt{ - (\lambda -a_1)(\lambda -a_2)(\lambda -a_3)(\lambda -a_4)} }}} ,\nonumber \\ I_{2,0} - i I_{2,1} - \frac{1}{2}{I_{1,0}}= & {} - \frac{1}{2}\int \limits _{{a_1}}^{{a_2}} {\frac{{\mathrm{d}\lambda }}{{{{(\lambda - i)}^2}\sqrt{ - (\lambda -a_1)(\lambda -a_2)(\lambda -a_3)(\lambda -a_4)} }}}. \end{aligned}$$

(44)

Then simple calculations render us

$$\begin{aligned} g_{1,1}= & {} {\text {Re}} C_1,\quad c_{1,1} = {\text {Im}} C_1,\quad C_1 = \frac{{{A_0}{C_0}{\alpha ^2}}}{{\alpha _1^2}},\nonumber \\ g_{1,3}= & {} \frac{{{A_0}{C_0}}}{2}\left( {1 - \frac{{{\alpha ^2}}}{{\alpha _1^2}}} \right) ,\quad c_{1,3} = - i g_{1,3},\nonumber \\ c_{2,1}= & {} - {\text {Re}}{C_2},\quad g_{2,1} = {\text {Im}}{C_2},\quad C_2 = \frac{{{A_0}C_0^2}}{{2\alpha _1^4}}\left[ {{\alpha ^4} + \frac{{{{(\alpha _1^2 - {\alpha ^2})}^2}}}{{2(\alpha _1^2 - 1)}}} \right] ,\nonumber \\ c_{2,2}= & {} - {\text {Re}}{C_3},\quad g_{2,2} = {\text {Im}}{C_3},\quad C_3 = \frac{{{A_0}C_0^2{{(\alpha _1^2 - {\alpha ^2})}^2}}}{{4\alpha _1^2(\alpha _1^2 - 1)(m - \alpha _1^2)}},\nonumber \\ c_{2,3}= & {} - \frac{{{A_0}C_0^2(\alpha _1^2 - {\alpha ^2})}}{{4\alpha _1^4}}\left[ {2{\alpha ^2} + \frac{{(2\alpha _1^2m + 2\alpha _1^2 - \alpha _1^4 - 3m)(\alpha _1^2 - {\alpha ^2})}}{{2(\alpha _1^2 - 1)(m - \alpha _1^2)}}} \right] ,\nonumber \\ g_{2,3}= & {} i{c_{2,3}}. \end{aligned}$$

(45)

With aid of equation (253.00) from Byrd and Friedman (1954), we also obtain:

$$\begin{aligned} c_{0,0} = A_0. \end{aligned}$$

(46)

1.1.2 Instance B

For integrals over the interval $(a_2,a_3)$, the only difference is in the parameters ${\alpha ^2}$, $\alpha _1^2$ and the parameter m. Thus, in (45) the values

$$\begin{aligned} \alpha ^2= \frac{a_3 - a_2}{a_3-a_1},\quad \alpha _1^2 = \frac{(a_3-a_2)(i-a_1)}{(a_3-a_1)(i-a_2)},\quad m = \frac{(a_3 - a_2)(a_4-a_1)}{(a_4-a_2)(a_3-a_1)} \end{aligned}$$

should be used.

1.1.3 Instance C

For integrals over the interval $(a_3,a_4)$ in (45), we have:

$$\begin{aligned} C_0= \frac{1}{a_3-i},\quad \alpha ^2= \frac{a_4-a_3}{a_4-a_2},\quad \alpha _1^2 = \frac{(a_4-a_3)(i-a_2)}{(a_4-a_2)(i-a_3)}. \end{aligned}$$

The parameters m and $A_0$ are the same as for the interval $(a_1,a_2)$.

Note The equality of m and $c_{0,0}$ in instances A and C means that, when the potential (4) has two minima, the periods of librations T about these on the same energy level are equal. This also holds true for the local minima corresponding to instances B and D, as the substitution $\lambda =\tan [(\varphi -{{\tilde{\varphi }}})/2]$ always allows us to get rid of the semi-infinite intervals of integration. Moreover, the values of $\left\langle \cos \varphi \right\rangle$ defined in (22) are also the same for librations around the two local minima. The values of $\left\langle \sin \varphi \right\rangle$ for librations about the two local minima differ by $\pi$, while the values of the adiabatic invariant (26) differ by $y/\sqrt{2}$.

1.1.4 Instance D

For integrals over two semi-infinite intervals, the values of the parameters in (45) are

$$\begin{aligned} C_0=\frac{1}{a_4-i},\quad \alpha ^2 = \frac{a_4-a_1}{a_3-a_1},\quad \alpha _1^2 = \frac{(a_4-a_1)(i-a_3)}{(a_3-a_1)(i-a_4)}. \end{aligned}$$

The parameter m is the same as for integration over $({a_2},{a_3})$, while $A_0$ is the same as for $({a_1},{a_2})$.

1.2 The case $a_1, a_2 \in {{\mathbb {R}}^1}$ ($a_1<a_2$), $a_3, a_4 \in {{\mathbb {C}}^1}$ ($a_3= {\bar{a}}_4$)

On this occasion, there are only two instances to be considered:

• A. Integration over $(a_1,a_2)$;

• B. Integration over $(-\infty ,a_1)\bigcup (a_2,+\infty )$.

1.2.1 Instance A

Here the following auxiliary parameters in the expressions for $c_{k,l}$ and $g_{k,l}$ will be used:

$$\begin{aligned}&\displaystyle \alpha = \frac{(a_1 A_2-a_2 A_1)-(A_2-A_1)i}{(a_1 A_2+a_2 A_1)-(A_2+A_1)i},\quad \alpha _1= \frac{A_2-A_1}{A_2+A_1},\\&\displaystyle C_0 = \frac{A_1+A_2}{(A_2 a_1-A_1 a_2)-(A_2-A_1)i}, \end{aligned}$$

where

$$\begin{aligned} {A_1}= & {} \sqrt{(a_1-a_0)^2 + b_0^2},\quad {A_2} = \sqrt{(a_2-a_0)^2 + b_0^2},\quad {A_0} = 1/\sqrt{A_1 A_2},\\ {a_0}= & {} {\text {Re}} a_3 = {\text {Re}} a_4,\quad {b_0} = {\text {Im}} a_3 = -{\text {Im}} a_4. \end{aligned}$$

Combining our Eq. (44) with formulae (259.04)⁴ and (341.01–341.04) from Byrd and Friedman (1954), we arrive at

$$\begin{aligned} c_{0,0}= & {} 2{A_0},\quad g_{1,1} = {\text {Re}} C_1,\quad c_{1,1} = {\text {Im}} C_1,\quad C_1 = 2 A_0 C_0 \alpha _1,\nonumber \\ c_{1,3}= & {} - \frac{i(\alpha - {\alpha _1})}{{(1 - {\alpha ^2})}}{A_0}{C_0},\quad g_{1,3} = i c_{1,3},\nonumber \\ c_{2,1}= & {} - {\text {Re}} C_2,\quad g_{2,1} = {\text {Im}} C_2,\quad C_2 = {A_0}C_0^2\left[ {\alpha _1^2 + \frac{{{{(\alpha - {\alpha _1})}^2}}}{{{\alpha ^2} - 1}}} \right] ,\nonumber \\ c_{2,2}= & {} -{\text {Re}} C_3,\quad g_{2,2} = {\text {Im}} C_3,\quad C_3 = \frac{{A_0 C_0^2 \alpha ^2{{(\alpha -\alpha _1)}^2}}}{{(1-\alpha ^2)(m + {\alpha ^2}{m'})}},\nonumber \\ c_{2,3}= & {} \frac{{{A_0}C_0^2}}{2}\frac{{(\alpha - {\alpha _1})}}{{({\alpha ^2} - 1)}}\left[ {2{\alpha _1} + \frac{{(\alpha - {\alpha _1})\left[ {{\alpha ^2}(2m - 1) - 2m} \right] }}{{({\alpha ^2} - 1)(m + {\alpha ^2}m')}}} \right] ,\quad {g_{2,3}} = i{c_{2,3}}. \end{aligned}$$

(47)

The parameter and characteristic of the elliptic integrals are calculated as follows:

$$\begin{aligned} m = \frac{{{{({a_2} - {a_1})}^2} - {{({A_2} - {A_1})}^2}}}{4 A_1 A_2},\quad n = \frac{{{\alpha ^2}}}{{{\alpha ^2} - 1}}. \end{aligned}$$

Here and hereafter, $m'$ denotes complimentary parameter $m'=1-m$.

1.2.2 Instance B

Expressions (47) stay the same. The auxiliary parameters are

$$\begin{aligned} \alpha= & {} \frac{{({a_2}{A_1} + {a_1}{A_2}) - ({A_2} + {A_1})i}}{{({a_2}{A_1} - {a_1}{A_2}) + ({A_2} - {A_1})i}},\quad {\alpha _1} = \frac{{{A_2} + {A_1}}}{{{A_1} - {A_2}}}, \\ C_{0}= & {} \frac{{{A_1} - {A_2}}}{{({A_2}{a_1} + {A_1}{a_2}) - ({A_2} + {A_1})i}}. \end{aligned}$$

The parameter m of the elliptic integrals is now given by

$$\begin{aligned} m = \frac{{{{({A_2} + {A_1})}^2} - {{({a_2} - {a_1})}^2}}}{{4{A_1}{A_2}}}. \end{aligned}$$

1.3 The case of ${a_j} \in {{\mathbb {C}}^1}$ $(j =\overline{1,4})$

Let us denote the real parts of the roots $a_j$ as $p_{1,2}$, and the imaginary parts as $b_{1,2}$:

$$\begin{aligned} a_{1,2}=p_1\pm b_1,\quad a_{3,4}=p_2\pm b_2. \end{aligned}$$

Without loss of generality, let us assume that ${p_1} < {p_2}$, ${b_1} > 0$, and ${b_2} > 0$ (if ${p_1} = {p_2}$, then also ${b_1} < {b_2}$). In the expressions for ${c_{k,l}}$ and ${g_{k,l}}$, the following auxiliary parameters will be employed:

$$\begin{aligned}&\displaystyle {A_1} = \sqrt{{{({p_2} - {p_1})}^2} + {{({b_2} - {b_1})}^2}},\quad {A_2} = \sqrt{{{({p_2} - {p_1})}^2} + {{({b_2} + {b_1})}^2}},\quad {A_0} = 2/({A_1} + {A_2}), \\&\displaystyle {g_1} = \sqrt{\frac{{4b_1^2 - {{({A_2} - {A_1})}^2}}}{{{{({A_2} + {A_1})}^2} - 4b_1^2}}},\quad \alpha = \frac{{{b_1} + {g_1}({p_1} - i)}}{{{p_1} - {b_1}{g_1} - i}},\quad {C_0} = \frac{1}{{{b_1} + {g_1}({p_1} - i)}}. \end{aligned}$$

Similarly to the previous case, we find:

$$\begin{aligned} c_{0,0}= & {} 2 A_0,\quad c_{1,1} = {\text {Im}} C_1,\quad g_{1,1} = {\text {Re}} C_1,\quad C_1 = \frac{2\alpha (1 + g_1 \alpha )}{1 + \alpha ^2} A_0 C_0,\\ g_{1,3}= & {} \frac{\alpha ^2(\alpha -g_1)A_0 C_0}{1+\alpha ^2},\quad c_{1,3} = -i g_{1,3},\quad c_{2,1} = -{\text {Re}} C_2,\quad g_{2,1} = {\text {Im}} C_2,\\ C_2= & {} A_0 C_0^2 \left[ g_1^2 + \frac{{2{g_1}(\alpha - {g_1})}}{{1 + {\alpha ^2}}} + \frac{{{{(\alpha - {g_1})}^2}}}{{(1 + {\alpha ^2})^2}} \right] ,\\ c_{2,2}= & {} {\text {Re}} C_3,\quad g_{2,2} = - {\text {Im}} C_3,\quad C_3 = \frac{{{A_0}C_0^2{{(\alpha - {g_1})}^2}{\alpha ^2}}}{{({\alpha ^2} + 1)({\alpha ^2} + m')}},\\ c_{2,3}= & {} -\frac{{{A_0}C_0^2}}{2}\frac{{(\alpha - {g_1}){\alpha ^2}}}{{({\alpha ^2} + 1)}}\left[ {2{g_1} + \frac{{(\alpha - {g_1})(2m' + 2{\alpha ^2} - {\alpha ^2}m)}}{{({\alpha ^2} + 1)({\alpha ^2} + m')}}} \right] ,\quad g_{2,3}= i c_{2,3}. \end{aligned}$$

The parameter and characteristic of the elliptic integrals are:

$$\begin{aligned} m = \frac{{4 {A_1}{A_2}}}{\left( {A_1} + {A_2}\right) ^2},\quad n = {\alpha ^2} + 1. \end{aligned}$$