On the coplanar eccentric non-restricted co-orbital dynamics

Abstract

We study the phase space of eccentric coplanar co-orbitals in the non-restricted case. Departing from the quasi-circular case, we describe the evolution of the phase space as the eccentricities increase. We find that over a given value of the eccentricity, around 0.5 for equal mass co-orbitals, important topological changes occur in the phase space. These changes lead to the emergence of new co-orbital configurations and open a continuous path between the previously distinct trojan domains near the $L_4$ and $L_5$ eccentric Lagrangian equilibria. These topological changes are shown to be linked with the reconnection of families of quasi-periodic orbits of non-maximal dimension.

Appendices

Appendix A: Time scales

The planar 3-body co-orbital problem has 3 time scales: the fast time scale, associated with the mean mean-motion $\eta ={\mathcal {O}} (1)$, the semi-fast time scale of fundamental frequency $\nu ={\mathcal {O}} (\sqrt{{\varepsilon } })$, associated with the evolution of the resonant angle $\zeta$, and the secular time scale of fundamental frequency $g_1$ and $g_2$, of order ${\mathcal {O}} ({\varepsilon } )$, associated with the evolution of the eccentricities and the arguments of perihelia. The separation of these time scales is a classical approach for the study of mean motion resonances (Henrard and Caranicolas 1989; Morbidelli 2002; Batygin and Morbidelli 2013; Delisle et al. 2012, 2014).

In theory, this separation allows for two averaging of the Hamiltonian: a first averaging over the fast angle $\lambda _2$ that we already considered in Sect. 2.3, and a second one over the semi-fast angle $\zeta$, in order to obtain the secular Hamiltonian. In the double averaged reduced case, we would obtain a 1 degree of freedom Hamiltonian which would describe the secular dynamics of the resonance. It would add an additional parameter $J_0$ (the action variable associated with the degree of freedom Z,$\zeta$). The canonical transformation for the variables $\varPi$ and ${\varDelta \varpi }$ associated with this second averaging differs from the identity only with coefficients of the order of ${\mathcal {O}} (\sqrt{\varepsilon } )$ (Morbidelli 2002).

1.1 A.1 Adiabatic invariants

In practice, this second averaging is rather difficult because the variables Z and $\zeta$ are not close to action-angle variables (Morais 1999, 2001; Beaugé and Roig 2001; Páez and Efthymiopoulos 2015). However, the possibility to do it gives us important information on the dynamics of the system: in the averaged reduced problem (2 degrees of freedom), the evolution of the variables $\varPi$ and ${\varDelta \varpi }$ is of size $\sqrt{{\varepsilon } }$ over durations of the order of $1/\nu$, these variables can be considered as constant on a time scale short with respect to 1 / g. For sufficiently low-mass co-orbitals, we can hence consider that the variables $\varPi$ and ${\varDelta \varpi }$ are adiabatic invariants.

1.2 A.2 Interpretation of numerical simulations

The change of coordinate (Eq. 15) from the variables of the planar 3-body problem to the variables of the averaged problem is ${\varepsilon }$ close from identity for all variables except $\zeta _2$. Similarly, the perturbations of the semi-fast time scale on the secular variables are of size $\sqrt{{\varepsilon } }$ (Morbidelli 2002). Thus, if we integrate numerically the full 3-body problem for co-orbital with low enough masses, for quasi-periodic orbits we can consider, on the one hand, the evolution of the variables ($Z,\zeta$) as their evolution in the averaged problem (they are ${\varepsilon }$ close), and on the other hand the evolution of the variables $e_j$ and $\varpi _j$ as their evolution in the secular problem (they are $\sqrt{\varepsilon }$ close).

Appendix B: Reference manifold

In this section, we aim to verify that all the trajectories of the phase space pass as close as we want from the reference manifold ${\mathcal {V}}$ defined by Eq. (35).

1.1 B.1 At first order in $e_j$

Equation (25) holds at first order in $e_j$. Hence, for $m_1=m_2$, all trajectories go through the plane $a_1=a_2$ twice per period $2\pi /\nu$. On the other hand, near a solution of the circular coplanar case $\zeta (t)$, the equation of variation in the direction ($x_j,{\tilde{x}}_j$), where the $x_j$ are the canonical Poincaré variables defined in Eq. (4), is given by the matrix (Robutel and Pousse 2013):

$$\begin{aligned} X= \begin{pmatrix} x_1\\ x_2 \end{pmatrix} \text {et} M(t)= i{\varepsilon } \eta \frac{m'_1 m'_2}{m_0} \begin{pmatrix} \frac{A(\zeta (t))}{m'_1} &{} \frac{\bar{B}(\zeta (t))}{\sqrt{m'_1 m'_2}}\\ \frac{B(\zeta (t))}{\sqrt{m'_1 m'_2}} &{} \frac{A(\zeta (t))}{m'_2} \end{pmatrix}, \end{aligned}$$

(36)

where A and B depend on the considered trajectory and on the time. For a given trajectory, since $\nu \gg g$, we can obtain an approximation of the secular dynamics in the direction ($x_j,{\tilde{x}}_j$) by averaging the expression of this matrix over a period $2\pi /\nu$ with respect to the time t. For equal mass co-orbitals, the symmetries of this matrix give relations of the form:

$$\begin{aligned} \begin{aligned} x_1&=\alpha \sqrt{m_2} {{\mathrm{e}}}^{i \left( \frac{\pi }{3}+g t\right) } + \beta \sqrt{m_1} {{\mathrm{e}}}^{i\frac{\pi }{3}}, \\ x_2&=- \alpha \sqrt{m_1} {{\mathrm{e}}}^{i g t} + \beta \sqrt{m_2}, \end{aligned} \end{aligned}$$

(37)

with $\alpha$ and $\beta$ complexes. Replacing these expression in the one of $\varPi$ (Eq. 32), and noting $\alpha \bar{\beta }= C {\text {e}}^{ic}$, we obtain:

$$\begin{aligned} \begin{aligned} \varPi&= (\alpha \bar{\alpha }-\beta \bar{\beta })(m_1-m_2)-2C\sqrt{m_1 m_2} \cos (g t + c). \end{aligned} \end{aligned}$$

(38)

On the other hand, we have:

$$\begin{aligned} \begin{aligned} {\varDelta \varpi }&=\arg (x_1\bar{x}_2), \ \ \text {where}\ \\ x_1\bar{x}_2&=\left[ \sqrt{m_1m_2}(\beta \bar{\beta }-\alpha \bar{\alpha })+Cm_2{{\mathrm{e}}}^{i(gt+c)}-Cm_1{{\mathrm{e}}}^{-i(gt+c)}\right] {{\mathrm{e}}}^{i\pi /3}. \end{aligned} \end{aligned}$$

(39)

When $m_1=m_2$, $\varPi$ librates around 0 with a frequency g. Using once more the expression (32), we obtain that the quantity $e_1^2-e_2^2$ behaves like an harmonic oscillator, librating around 0 with the frequency $g=2|c|={\mathcal {O}} ({\varepsilon } )$. All the trajectories of the phase space hence go through the plane $e_1=e_2$ twice per period $2\pi /g$.

As long as $\nu$ and g are non-resonant, all trajectories get as close as we want to the manifold ${\mathcal {V}}$ in a finite time.

Fig. 17

Minimal value of the quantity $\log \left( (a_1/\bar{a}-a_2/\bar{a})^2+(e_1-e_2)^2 \right)$ over $10\times 10^{5}$ orbital periods with a step of 0.01 orbital period and a fixed value of the angular momentum $J_1(e_1=e_2=0.4)$, with the following initial conditions: $a_1=a_2=1$, $m_1=m_2=10^{-5} m_0$ for the top row and $m_2=3m_1=1.5\, 10^{-5}$ for the bottom one. a $e_1=0.00$ and $e_2 \approx 0.55$; b $e_1=0.10$ and $e_2 \approx 0.54$; c $e_1=0.20$ and $e_2 \approx 0.52$; d $e_1=0.30$ and $e_2 \approx 0.47$. The trajectories ejected before the end of the integration are identified by white pixels. See Sect. 4 for more details about the integrations

1.2 B.2 Large eccentricities

We check numerically if the definition ${\mathcal {V}} =\{a_1=a_2,e_1=e_2\}$ holds for higher eccentricities. Note that we consider only trajectories that reach $a_1=a_2$ on their orbit. To perform this check, we take grids of initial conditions for $\zeta \in [ 0^\circ : 360^\circ ]$ and ${\varDelta \varpi } \in [-180^\circ :180^\circ ]$ and several values of $\varPi$ for a fixed value of $J_1$ such that $J_1=J_1(e_1=e_2=0.4)$. The corresponding values of the eccentricities are given by:

$$\begin{aligned} \begin{aligned} e_2=\sqrt{1-\frac{2J_1}{\varLambda ^0_1}-\sqrt{1-e_1^2}}. \end{aligned} \end{aligned}$$

(40)

Figure 17 shows the value of $(a_1/\bar{a}-a_2/\bar{a})^2+(e_1-e_2)^2$ for several values of $\varPi$ for $m_1=m_2$ (top line) and $m_1 \ne m_2$ (bottom line). The integrations are conducted over $10/{\varepsilon }$ orbital periods, hence only a few times $2\pi /g$ at best. For all initial conditions when $m_1=m_2$, the criterion

$$\begin{aligned} \frac{(a_1-a_2)^2}{\bar{a}^2}+(e_1-e_2)^2 < \epsilon _\varSigma \end{aligned}$$

(41)

is met for $\epsilon _\varSigma \approx 10^{-8}$. Although this verification is not exhaustive, it suggests that the chosen reference manifold represents a significant part of the phase space of the averaged reduced problem. However it is possible that, especially at high eccentricities, stable domains appear for which the orbits never reach $e_1=e_2$ even in the case $m_1=m_2$, but none was discovered during this study. A study performed in the case $e_1=e_2=0.7$ yielded similar results.

In order to compare with the case $m_1 \ne m_2$, the bottom line of Fig. 17 shows that there are areas of the phase space where the criterion (41) is not verified for $\epsilon '_\varSigma =10^{-4}$. There are hence orbits in the phase space that are not represented by the trajectories taking their initial conditions on the manifold ${\mathcal {V}} =\{a_1=a_2,e_1=e_2\}$. This is not surprising: we know, for example, that, at least for moderate eccentricities, the position of the $AL_4$ equilibrium is approximated by $m_1e_1=m_2e_2$. In the case $m_1 \ne m_2$, any trajectory librating sufficiently close to this equilibrium would never cross the $e_1=e_2$ manifold.

Appendix C: Identification of the ${\mathcal {F}}$ families

We show here how the separation of the time scales allows us to identify the position of the ${\mathcal {F}}$ anywhere in the phase space.

1.1 C.1 Identification of the ${{\mathcal {F}} ^{sc}}$ families

The ${{\mathcal {F}} ^{sc}}$ families are families of periodic orbit of the reduced averaged problem, whose period is associated with the secular time scale. The position of the ${{\mathcal {F}} ^{sc}}$ families can be identified by studying the critical points of the averaged Hamiltonian. Let us use the hypothesis of adiabatic invariant for the variables $\varPi$ and ${\varDelta \varpi }$ (see Appendix A): on a short time scale with respect to 1 / g, ${{\mathcal {F}} ^{sc}}$ is made of orbits that behave as fixed points of the reduced averaged problem. The orbits belonging to ${{\mathcal {F}} ^{sc}}$ are thus orbits which satisfy:

$$\begin{aligned} \frac{\partial }{\partial Z}{\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} } = \frac{\partial }{\partial \zeta } {\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }= 0. \end{aligned}$$

(42)

where Z and $\zeta$ are conjugated canonical variables. This is equivalent to:

$$\begin{aligned} \dot{\zeta }= \dot{Z}= 0. \end{aligned}$$

(43)

Starting from the reduced Hamiltonian (Sect. 3.1), we can estimate the value of the averaged Hamiltonian at any point of the phase space by doing a numerical averaging over the fast angle Q. We can identify the orbits belonging to ${{\mathcal {F}} ^{sc}}$ by finding the orbits for which $\dot{Z} = 0$ on the manifold $\dot{\zeta }=0$. For a given value of the constants $\varPi$ and ${\varDelta \varpi }$, we take a grid of values for $\zeta$ and estimate the averaged Hamiltonian at each point. We can then have the approximate position of the points where $\dot{Z} = 0$ by finding the positions on the grid where the equation

$$\begin{aligned} \frac{\partial }{\partial \zeta } {\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }|_{\zeta =\zeta _k} \times \frac{\partial }{\partial \zeta } {\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }|_{\zeta =\zeta _{k+1}} < 0 \end{aligned}$$

(44)

is satisfied, with

$$\begin{aligned} \frac{\partial }{\partial \zeta } {\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }|_{\zeta =\zeta _k} = \frac{{\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }(Z,\zeta _{k+1},{\varDelta \varpi } ,\varPi ) -{\mathcal {H}} _{{\mathcal {R}} {\mathcal {M}} }(Z,\zeta _{k-1},{\varDelta \varpi } ,\varPi ) }{|\zeta _{k+1}-\zeta _{k-1}|}. \end{aligned}$$

(45)

Note that it is not guaranteed that the associated trajectory in the full 3-body problem is quasi-periodic.

Alternatively, numerical integrations allow us to determine an empiric criterion for a numerical determination of the position of ${{\mathcal {F}} ^{sc}}$. In the various integrations that we computed through this study, we noted that the amplitude of variation of Z (hence $a_1-a_2$) seems not to be impacted much by the frequency g. We hence make the hypothesis that if an orbit in a regular area of the phase space verifies the condition

$$\begin{aligned} (\max {(Z)}-\min {(Z)} ) < \epsilon _\nu , \end{aligned}$$

(46)

with $\epsilon _\nu \propto \sqrt{{\varepsilon } }$, this orbit is in the neighbourhood of the manifold ${{\mathcal {F}} ^{sc}}$. One can check in Figs. 5, 6, 7, 8, 9, 10, 11, 12, 13 and 14, where the quasi-periodic orbits that verify Eq. (46) are identified by brown pixels, and the point of the phase space satisfying Eq. (42) is identified by purple dots that both methods yield very similar results in the regular area of the phase space.

1.2 C.2 Application in the case $m_1=m_2$

We can apply the research of the critical points of the Hamiltonian to identify the position of ${{\mathcal {F}} ^{sc}}$ in the case $m_1=m_2$. We assume that, as it is the case for circular co-orbitals, the manifold $\dot{\zeta }=0$ is located at $Z=0$. We know that the equilibriums $L_k$ and $AL_k$ are all located in the plane $\varPi =0$ ($e_1=e_2$). We can hence explore the manifold ${\mathcal {V}} =\{Z,\zeta , \varPi ,{\varDelta \varpi } / Z=\varPi =0 \}$. We chose a grid of initial condition for $\zeta$ and ${\varDelta \varpi }$ with a step of $0.5^\circ$, and we compute numerically the averaged Hamiltonian at each point of the grid. In Fig. 4, we show all the points of ${\mathcal {V}}$ that verify the condition (44). Each graph corresponds to a different value of the total angular momentum (different value of $e_1=e_2$).

1.3 C.3 Identification of the ${{\mathcal {F}} ^{{ sf}}}$ families

The ${{\mathcal {F}} ^{{ sf}}}$ families are families of periodic orbits of the reduced averaged problem, whose period is associated with the semi-fast (resonant) time scale. The method developed in Sect. C.1 cannot be used directly to determine the position of the ${{\mathcal {F}} ^{{ sf}}}$ manifold because it requires to numerically average the Hamiltonian over the semi-fast angle $\zeta$, which is somehow laborious, see Sect. A.

However, the evolution of the variables ${\varDelta \varpi }$ and $\varPi$ during the numerical integrations of the 3-body problem is ${\mathcal {O}} (\sqrt{\varepsilon } )$ close to their evolution in the secular problem (see Sect. A). Since the orbits belonging to ${{\mathcal {F}} ^{{ sf}}}$ are fixed points of the 1-degree of freedom secular problem, we make the following hypothesis: all orbits in a regular area of the phase space (far from the separatrix, the chaotic and the unstable areas) that verify

$$\begin{aligned} (\max {({\varDelta \varpi } )}-\min {({\varDelta \varpi } )} ) < \epsilon _g, \end{aligned}$$

(47)

with $\epsilon _g \propto \sqrt{{\varepsilon } }$ are in the neighbourhood of ${{\mathcal {F}} ^{{ sf}}}$. One can check that such orbits (represented by black pixel in Figs. 5, 6, 7, 8, 9, 10, 11, 12, 13 and 14) are indeed in the neighbourhood of the analytical approximation of the positions of the ${{\mathcal {F}} ^{{ sf}}}$ families, see Leleu (2016, Sect. 2.7.2).