© ESO 2018

1 Introduction

In studying the dynamics of satellites in the solar system, it is important to establish whether the bodies are bound to the planets or may escape to heliocentric orbits. Therefore, the perturbation by the Sun is a key issue when the stability of these satellites is to be addressed. In the case of the Moon, the development of Hill’s lunar theory allowed identifying the so-called evection term (Brouwer & Clemence 1961) as the largest periodic correction to the mean longitude of the Moon. For any satellite, the evection term is associated with the harmonic cos(2λ_⊙− 2ϖ), where λ_⊙ is the mean longitude of the Sun and ϖ is the longitude of the pericenter of the satellite. The evection resonance arises when the precession rate of ϖ equals the solar mean motion, causing the angle λ_⊙−ϖ to librate around an equilibrium point. In the classical evection, the pericenter precession is driven by the solar perturbation itself, but in principle, any other perturbation inducing a pericenter precession may originate an evection resonance.

In the restricted three-body problem, Henon (1969, 1970) showed that the classical evection resonance appears as a bifurcation of a family of simple periodic orbits at a value of the semi-major axis a = 0.45 R_Hill.

Hamilton & Krivov (1997) showed that in the case of prograde orbits, the evection resonance appears at a = 0.53 R_Hill and is characterized by the resonant angle λ_⊙−ϖ librating either around 0° or 180°. This alignment or anti-alignment of the satellite pericenter with the Sun direction induces cumulative perturbations that may cause the escape of the satellite (Nesvorný et al. 2003), which might be relevant for the stability of bodies that are migrating as a result of tidal evolution or gas drag.

The evection resonance plays an important role in sculpting the architecture of the satellite systems. Nesvorný et al. (2003) studied the orbital and collisional evolution of satellites and found that prograde satellite orbits with large semi-major axes are unstable because of the effect of the evection resonance. Ćuk & Gladman (2009) explored the fate of fictitious objects trapped in the lunar Trojan points after the formation of the Moon, and found that these bodies can survive the tidal migration of the Moon until they reach 38 Earth radii, where the evection resonance ejects them from the system. Ćuk et al. (2016) also studied the past evolution of the Tethys–Dione system and found that the evection resonance perturbing a pair of medium-size moons is the most likely mechanism for triggering instability and massive collisions in the Saturn system. Finally, Spalding et al. (2016) studied the inward migration of giant exoplanets immersed in protoplanetary gas disks, and showed that dynamical interactions may naturally destroy the hypothetical moons of such planets through capture into the evection resonance. Several analytical approaches to the study of the evection resonance have been published. Yokoyama et al. (2008) analytically determined that the semi-major axis for the appearance of the evection resonance lies at a = 0.529R_Hill, while Frouard et al. (2010) used dynamical maps to show that this resonance constitutes the outermost region of stability for the prograde satellites of Jupiter, being surrounded by chaotic orbits. The authors concluded that the evection can harbor stable orbits at distances from Jupiter of ~0.43 au. Frouard et al. (2010) also showed that the inclusion of third-order terms in the expansion of the solar disturbing function in Legendre polynomial causes the evection resonance to display a non-symmetric topology with respect to λ _⊙−ϖ = 0° and 180°. Andrade-Ines & Robutel (2018) constructed a secular theory for planetary motion around compact binaries and analyzed the occurrence of the evection resonance that is driven by the stellar companion. The study of the evection resonance using analytical models is not straightforward, however, because usually large perturbations act on the regions of interest. In this work, we study the occurrence and stability of the evection resonance in the satellite system of Saturn. This system is particularly interesting because it harbors the only known examples of co-orbital configurations in satellites in the solar system. These are

• 1.

  the pair Janus and Epimetheus, which evolve in a relative horseshoe orbit (Murray & Dermott 1999), and may have originated from a former Trojan configuration (Gott 2005);

• 2.

  the Trojan system of Tethys, accompanied by the moons Telesto (at L₄) and Calypso (at L₅) (see Niederman et al. 2018, and references therein);

• 3.

  the Trojan system of Dione, accompanied by the moons Helene (at L₄) and Polydeuces (at L₅).

Our main goal is to provide some clues about the possible past dynamics of both single and co-orbital moons. In Sect. 2 we introduce an analytical expansion of the Hamiltonian for two satellites in Trojan co-orbital configuration, and we apply it to describe the basic topology of the evection resonance in such a system. In Sect. 3 we introduce our numerical analysis of the stability of the evection resonance in the Saturn system, and we apply it first to the case of a single satellite (Sect. 3.1), and then to the case of a pair of Trojan moons (Sect. 3.2). Section 4 is devoted to studying the effect of the evection resonance on satellites that migrate as a result of tidal evolution. Finally, our conclusions are presented in Sect. 5.

2 Analytical model

In this section, we set up a model to describe the evection resonance in the case of two satellites in co-orbital Trojan motion. We then apply this model to the Saturn system.

2.1 Trojan pair under the influence of a distant perturber

We consider a hierarchical system of four bodies of masses m_i (i = 0, ..., 3), in which m₀ is a central oblate body of radius R₀ (planet), m₁, m₂ ≪ m₀ are two co-orbital bodies around m₀ (satellites), and m₃ is a distant perturber (the Sun). u_i, w_i are the barycentric positions and momenta of the bodies, and we introduce the following set of canonical coordinates: $$\begin{array}{l}{\text{r}}_0=\frac{1}{M}{\sum }_{j=0}^3{m}_j{\text{u}}_j,\\ {\text{r}}_i={\text{u}}_i-{\text{u}}_0,i=1,2,\\ {\text{r}}_3={\text{u}}_3-\frac{1}{m}{\sum }_{j=0}^2{m}_j{\text{u}}_j,\end{array}$$(1)

where $m={\sum }_{i=0}^2{m}_i$ and M = m + m₃. The positions r ₁, r ₂ are measured from m₀, while r ₃ is measured from the barycenter of the sub-system m₀, m₁, m₂. Their conjugated momenta are $$\begin{array}{l}{\text{p}}_0={\sum }_{j=0}^3{\text{w}}_j,\\ {\text{p}}_i={\text{w}}_i-\frac{m_i}{m}{\sum }_{j=0}^2{\text{w}}_j,i=1,2,\\ {\text{p}}_3={\text{w}}_3-\frac{m_3}{M}{\sum }_{j=0}^3{\text{w}}_j.\end{array}$$(2)

Defining $$\begin{array}{lll}{\mu }_i=\frac{m_i{m}_0}{m_i+m_0},\hfill & \hfill & {\beta }_i=G(m_0+m_i),i=1,2,\hfill \\ {\mu }_3=\frac{m_{3}m}{m_3+m},\hfill & \hfill & {\beta }_3=GM,\hfill \end{array}$$(3)

the Hamiltonian of the system reads $$\mathcal{H}={\sum }_{i=1}^3\left(\frac{p_i^2}{2{\mu }_i}-\frac{{\beta }_i{\mu }_i}{r_i}\right)+{\sum }_{i=1}^2{\mathcal{J}}_i+{\mathcal{R}}_{12}+{\mathcal{R}}_3,$$(4)

where ${\mathcal{J}}_i$ is the oblateness perturbation felt by m_i, $${\mathcal{R}}_{12}=-\frac{G{m}_1{m}_2}{\left|{\text{r}}_1-{\text{r}}_2\right|}+\frac{1}{m_0}{\text{p}}_1\cdot {\text{p}}_2$$(5)

is the co-orbital disturbing function, and $$\begin{array}{lll}{\mathcal{R}}_3\hfill & =\hfill & G{m}_3{m}_0\left(\frac{1}{r_3}-\frac{1}{r_3^{\prime }}\right)\hfill \\ \hfill & \hfill & +G{m}_3{\sum }_{i=1}^2{m}_i\left(\frac{1}{r_3}-\frac{1}{\left|{\text{r}}_3^{\prime }-{\text{r}}_i\right|}\right)\hfill \end{array}$$(6)

is the distant perturbation, with the auxiliary variable $${\text{r}}_3^{\prime }={\text{r}}_3+\frac{1}{m}{\sum }_{i=1}^2{m}_i{\text{r}}_i,$$(7)

which is the position of m₃ relative to m₀ (Chambers 2010).

Assuming that R₀ ≪ r_i (i = 1, 2), we may expand ${\mathcal{J}}_i$ in powers of R₀∕r_i using Legendre polynomials. Up to the fourth order, we obtain $$\begin{array}{lll}{\mathcal{J}}_i\hfill & =\hfill & {\beta }_i{\mu }_i{J}_2\frac{R_0^2}{r_i^3}\left(\frac{3}{2}{\sin }^2{\phi }_i-\frac{1}{2}\right)\hfill \\ \hfill & \hfill & +{\beta }_i{\mu }_i{J}_4\frac{R_0^4}{r_i^5}\left(\frac{35}{8}{\sin }^4{\phi }_i-\frac{30}{8}{\sin }^2{\phi }_i+\frac{3}{8}\right)+\cdots \hfill \end{array}$$(8)

where J_k are the zonal harmonic coefficients of m₀, and φ_i is the latitude of m_i over the m₀ equator (Murray & Dermott 1999).

Similarly, assuming that r_i ≪r₃ (i = 1, 2), we may expand ${\mathcal{R}}_3$ in powers of r_i∕r₃. Taking into account Eq. (7) and introducing the small quantities $$\overline{\nu }=\frac{m_1{m}_2}{m},{\nu }_i=\frac{m_i}{m},i=1,2,$$(9)

we obtain up to the third order¹ $${\mathcal{R}}_3={\mathcal{R}}_3^{(2)}+{\mathcal{R}}_3^{(2c)}+{\mathcal{R}}_3^{(3)}+{\mathcal{R}}_3^{(3c)}+\cdots ,$$(10)

with $$\begin{array}{lll}\hfill & \hfill & {\mathcal{R}}_3^{(2)}=-G{m}_3{\sum }_{i=1}^2{m}_i\left(1-{\nu }_i\right)\frac{r_i^2}{r_3^3}\left(\frac{3}{2}{\cos }^2{\gamma }_i-\frac{1}{2}\right),\hfill \\ \hfill & \hfill & {\mathcal{R}}_3^{(2c)}=G{m}_3\overline{\nu }\frac{r_1{r}_2}{r_3^3}\left(3\cos {\gamma }_1\cos {\gamma }_2-\cos \varphi \right),\hfill \\ \hfill & \hfill & {\mathcal{R}}_3^{(3)}=-G{m}_3{\sum }_{i=1}^2{m}_i\left(1-3{\nu }_i+2{\nu }_i^2\right)\hfill \\ \hfill & \hfill & \times \frac{r_i^3}{r_3^4}\left(\frac{5}{2}{\cos }^3{\gamma }_i-\frac{3}{2}\cos {\gamma }_i\right),\hfill \\ \hfill & \hfill & {\mathcal{R}}_3^{(3c)}=G{m}_3\overline{\nu }\left(1-2{\nu }_1\right)\frac{r_1^2{r}_2}{r_3^4}(\frac{15}{2}{\cos }^2{\gamma }_1\cos {\gamma }_2\hfill \\ \hfill & \hfill & -\frac{3}{2}\cos {\gamma }_2-3\cos {\gamma }_1\cos \varphi )\hfill \\ \hfill & \hfill & +G{m}_3\overline{\nu }\left(1-2{\nu }_2\right)\frac{r_1{r}_2^2}{r_3^4}(\frac{15}{2}{\cos }^2{\gamma }_2\cos {\gamma }_1-\frac{3}{2}\cos {\gamma }_1\hfill \\ \hfill & \hfill & -3\cos {\gamma }_2\cos \varphi \text{ vphantom}\frac{15}{2}),\hfill \end{array}$$

where γ_i is the angle between r₃ and r _i and ϕ is the angle between r₁ and r ₂. The terms ${\mathcal{R}}_3^{(ic)}$ arise from the motion of the co-orbital pair. When either m₁ or m₂ is zero, then ${\mathcal{R}}_3^{(ic)}=0$ and Eq. (10) reduces to the classical expansion of the lunar theory (Brouwer & Clemence 1961).

No expansion is made of the co-orbital perturbation, since r₁ ~ r₂, but taking into account that the osculating velocity of m_i is $$\begin{array}{ll}{\text{v}}_i=\frac{{\text{p}}_i}{{\mu }_i}=\frac{m_0+m_i}{m_{0}m}\left[(m_0+m_j){\stackrel{.}{\text{r}}}_i-m_j{\stackrel{.}{\text{r}}}_j\right],\hfill & i=1,2,\hfill \\ \hfill & j\ne i,\hfill \end{array}$$(15)

we may write $$\begin{array}{lll}{\mathcal{R}}_{12}\hfill & =\hfill & -G{m}_1{m}_2{\left(r_1^2+r_2^2-2r_1{r}_2\cos \varphi \right)}^{-1/2}\hfill \\ \hfill & \hfill & +\overline{\nu }\left(1+\frac{2\overline{\nu }}{m_0}\right){\dot{r}}_1{\dot{r}}_2\cos {\varphi }^{\prime }-\frac{{\overline{\nu }}^2}{{\mu }_1{\mu }_2}\left({\mu }_1{\dot{r}}_1^2+{\mu }_2{\dot{r}}_2^2\right)\hfill \end{array}$$(16)

where ϕ′ is the angle between ${\stackrel{.}{\text{r}}}_1$ and ${\stackrel{.}{\text{r}}}_2$.

Equations (8), (10), and (16) can be further expanded in terms of the orbital elements of the Keplerian motion. We restrict the motion to coplanar orbits and let the equator of m₀ coincide with the reference plane (φ_i = 0). We introduce the set of canonical elements: $$\begin{array}{lll}\hfill & \hfill & {\theta }_1=\frac{{\lambda }_2-{\lambda }_1}{2};Z_1={\mu }_2\sqrt{{\beta }_2{a}_2}-{\mu }_1\sqrt{{\beta }_1{a}_1}\hfill \\ \hfill & \hfill & {\theta }_2=\frac{{\lambda }_2+{\lambda }_1}{2};Z_2={\mu }_2\sqrt{{\beta }_2{a}_2}+{\mu }_1\sqrt{{\beta }_1{a}_1}\hfill \\ \hfill & \hfill & {\theta }_3={\lambda }_3;Z_3={\mu }_3\sqrt{{\beta }_3{a}_3}-W_1-W_2\hfill \\ \hfill & \hfill & {\psi }_i={\lambda }_3-{\varpi }_i;W_i={\mu }_i\sqrt{{\beta }_i{a}_i}\left(1-\sqrt{1-e_i^2}\right),i=1,2,\hfill \\ \hfill & \hfill & {\psi }_3=-{\varpi }_3;W_3={\mu }_3\sqrt{{\beta }_3{a}_3}\left(1-\sqrt{1-e_3^2}\right),\hfill \end{array}$$(17)

where a_i, e_i, λ_i, and ϖ_i are the semi-major axis, eccentricity, mean longitude, and pericenter longitude, respectively, and we expand the Hamiltonian in powers of the eccentricities, around e_i = 0, up to $\mathcal{O}(e^4)$.

The co-orbital motion is characterized by the libration of the resonant angle θ ₁, while the synodic angle θ₂ is a fast angle that can be eliminated by an averaging: $$\overline{\mathcal{H}}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\mathcal{H}}\text{d}{\theta }_2,$$(18)

thus Z₂ becomes an integral of motion.

In this way, the averaged co-orbital perturbation, ${\overline{\mathcal{R}}}_{12}$, reduces to an expansion similar to that obtained by Robutel & Pousse (2013). This expansion has a singularity when a₁ = a₂ and θ ₁ = 0, even if e_i ≠ 0, and only contains terms of degree 1, e², e⁴, since the odd-degree terms are eliminated by the average. The averaged expansion of the oblateness perturbation, ${\overline{\mathcal{J}}}_i$, is straightforward and also contains terms of degree 1, e², e⁴. Only the averaged expansion of the distant perturbation, ${\overline{\mathcal{R}}}_3$, contains both even- and odd-degree terms in e. In particular, ${\overline{\mathcal{R}}}_3^{(2)}$ only contains terms of degree 1, e², while ${\overline{\mathcal{R}}}_3^{(2c)}$ contains terms of degree 1, e², e⁴. Both ${\overline{\mathcal{R}}}_3^{(3)}$ and ${\overline{\mathcal{R}}}_3^{(3c)}$ contain terms of degree e, e³ only.

After rearranging terms, the averaged Hamiltonian may be written as follows: $$\begin{array}{lll}\overline{\mathcal{H}}\hfill & =\hfill & {\overline{\mathcal{H}}}_0\left({\theta }_1,Z_1,Z_2,Z_3\right)\hfill \\ \hfill & \hfill & +{\overline{\mathcal{H}}}_1\left({\theta }_1,{\psi }_1,{\psi }_2,Z_1,Z_2,Z_3,W_1,W_2\right)\hfill \\ \hfill & \hfill & +{\overline{\mathcal{H}}}_2\left({\theta }_1,{\theta }_3,{\psi }_1,{\psi }_2,{\psi }_3,Z_1,Z_2,Z_3,W_1,W_2,W_3\right)\hfill \end{array}$$(19)

where ${\overline{\mathcal{H}}}_0$ is a one degree of freedom Hamiltonian that contains all the terms that do not depend on the eccentricities, ${\overline{\mathcal{H}}}_1$ is a perturbation that contains the terms that depend solely on e₁, e₂, and ${\overline{\mathcal{H}}}_2$ is a perturbation that contains all the termsdepending on e₃. Elimination of two degrees of freedom is immediate by setting e₃ = 0. In that case, ${\overline{\mathcal{H}}}_2$ vanishes and the resulting Hamiltoniandoes not depend on W₃, ψ₃, nor on θ₃, thus Z₃ becomes an integral of motion. This integral reflects the perturbation of the co-orbital system onto the distant perturber, which in this approach has a fixed circular orbit with a varying semi-major axis.

The literal expansion of $\overline{\mathcal{H}}$ in non-canonical elements contains more than 50 terms, but we can restrict it to those of the lowest degree and to the most relevant terms. For the distant perturbation, we restrict the terms arising from the co-orbital motion, ${\overline{\mathcal{R}}}_3^{(2c)}$ and ${\overline{\mathcal{R}}}_3^{(3c)}$, to degree e², and keep the main terms, ${\overline{\mathcal{R}}}_3^{(2)}$ and ${\overline{\mathcal{R}}}_3^{(3)}$, up to degree e³. This is a reasonable assumption when we take into account that ${\nu }_i\lesssim \mathcal{O}(r_i/r_3)$ for most satellite systems. Therefore, the term ${\overline{\mathcal{R}}}_3^{(3c)}$ is $\mathcal{O}({\nu }_i)$ smaller than the term ${\overline{\mathcal{R}}}_3^{(3)}$, while the term ${\overline{\mathcal{R}}}_3^{(2c)}$ is on the same order as ${\overline{\mathcal{R}}}_3^{(3)}$. For the co-orbital perturbation, ${\overline{\mathcal{R}}}_{12}$, we follow Robutel & Pousse (2013) and keep only terms of degree e². Finally, for the oblateness perturbation, we assume J₄ = 0. The averaged Hamiltonian then reduces to the following relevant terms: $$\begin{array}{lll}\overline{\mathcal{H}}\hfill & =\hfill & -{\sum }_{i=1}^3\frac{{\beta }_i{\mu }_i}{2a_i}+{\sum }_{i=1}^2\left(\frac{A_{0i}}{a_3^3}+D_i\right)+B_0+\frac{C_0}{a_3^3}\hfill \\ \hfill & \hfill & +{\sum }_{i=1}^{2}3\frac{A_{1i}}{a_3^4}{e}_i\cos {\psi }_i+{\sum }_{i=1}^2\frac{3}{2}\left(\frac{A_{0i}}{a_3^3}+D_i\right)e_i^2\hfill \\ \hfill & \hfill & +{\sum }_{i=1}^2\left(B_1-\frac{1}{2}\frac{C_0}{a_3^3}\right)e_i^2+{\sum }_{i=1}^2\frac{15}{8}{D}_i{e}_i^4\hfill \\ \hfill & \hfill & +{\sum }_{i=1}^2\frac{3}{2}\left(5\frac{A_{0i}}{a_3^3}+\frac{1}{4}\frac{C_0}{a_3^3}\right)e_i^2\cos 2{\psi }_i\hfill \\ \hfill & \hfill & +\left(B_2+\frac{C_1+C_2}{a_3^3}\right)e_1{e}_2\cos ({\psi }_1-{\psi }_2)\hfill \\ \hfill & \hfill & +3\frac{C_2}{a_3^3}{e}_1{e}_2\cos ({\psi }_1+{\psi }_2)\hfill \\ \hfill & \hfill & +\left(B_3+\frac{C_3}{a_3^3}\right)e_1{e}_2\sin ({\psi }_1-{\psi }_2)\hfill \\ \hfill & \hfill & +\frac{C_4}{a_3^3}\left(e_1^2\sin 2{\psi }_1-e_2^2\sin 2{\psi }_2\right)\hfill \\ \hfill & \hfill & +{\sum }_{i=1}^2\frac{1}{4}\frac{A_{1i}}{a_3^4}{e}_i^3\left(9\cos {\psi }_i+35\cos 3{\psi }_i\right)\hfill \end{array}$$(20)

where A_ki, B_k, C_k, and D_i are functions of θ₁, Z₁, and Z₂. Expressions for these functions in terms of non-canonical elements are given in the appendix.

We note that when m₃ = 0, all the A_ki and C_k functions are zero and theexpansion reduces to the same expansion as reported by Robutel & Pousse (2013), with the addition of the oblateness terms D_i. On the other hand, when m₂ = 0, all the B_k and C_k functions are zero, and the expansion reduces to the classical expansion of the lunar theory, again with the addition of the oblateness terms (Frouard et al. 2010).

For e₁ = e₂ = 0, the Hamiltonian $$\begin{array}{ll}{\overline{\mathcal{H}}}_0\hfill & =-{\sum }_{i=1}^3\frac{{\beta }_i{\mu }_i}{2a_i}+{\sum }_{i=1}^2\left(\frac{A_{0i}}{a_3^3}+D_i\right)+B_0+\frac{C_0}{a_3^3}\hfill \end{array}$$(21)

reflects the zeroth-order topology of the co-orbital motion. This is shown in Fig. 1 (left panel) for the Dione–Helene system. We note that at variance with the model of Robutel & Pousse (2013), the inclusion in ${\overline{\mathcal{H}}}_0$ of the zeroth-order terms of the solar perturbation (A_0i, C₀) and the planet oblateness (D_i) shift the libration center to values lower than 60° at increasing distances from the planet, as shown in Fig. 1 (right panel).

The Hamiltonian Eq. (20) represents a system with three degrees of freedom that might be reduced in principle, provided that the degrees of freedom evolve on very different timescales. This would allow the application of the adiabatic invariance principle to eliminate degrees of freedom, as in Beaugé & Roig (2001) or Saillenfest et al. (2016, 2017). This rigorous semi-analytical treatment is beyond the scope of this work. Here, we instead use a very simple zeroth-order approach to provide an overview of the Hamiltonian topology. We know that in the absence of other perturbations, the stable equilibrium of the Trojan co-orbital motion arises when Δ λ = λ₂ −λ₁ ≃ 60°, Δ ϖ= ϖ₂ −ϖ₁ ≃ 60°, and e₂ ≃ e₁ (Giuppone et al. 2010). This corresponds to the so-called asymmetric apsidal corotation resonance solution, or ACR, which is schematically represented in Fig. 2. Therefore, we may freeze two degrees of freedom in Eq. (20) by fixing θ ₁, Z₁ at one of the stable equilibrium points of ${\overline{\mathcal{H}}}_0$, and also setting ψ₁ −ψ₂ = 60° and e₁ = e₂. In this way, the functions A_ki, B_k, C_k, and D_i become constant coefficients. The resulting level curves of the remaining Hamiltonian of one degree of freedom are shown in Figs. 3–6.

It is worth recalling that in freezing the system at the ACR solution, we avoid the occurrence of other possible equilibria of the Hamiltonian with three degrees of freedom in Eq. (20). Moreover, we cannot be sure a priori that the ACR solution represents an equilibrium solution when the solar perturbation and the oblateness perturbations are taken into account. Nevertheless, the numerical experiments presented in Sects. 3.2 and 4 demonstrate that the ACR assumption is still valid in this case, and the ACR is an effective equilibrium state of the system. Exploring other equilibria would require searching for the zero-amplitude solutions given by the stationary conditions: $$\frac{\partial \overline{\mathcal{H}}}{\partial {\Theta }_i}=\frac{\partial \overline{\mathcal{H}}}{\partial J_i}=0,$$(22)

where J_i are the actions (Z₁, W₁, W₂) and Θ_i are the angles (θ₁, ψ₁, ψ₂). These solutions represent the extrema of Eq. (20) in a six-dimensional space and could be assessed using numerical (e.g., Giuppone & Leiva 2016), semi-analytical (e.g., Giuppone et al. 2010), or analytical (e.g., Robutel & Pousse 2013) techniques. Such an analysis, however, is beyond the scope of this work.

Fig. 1 Left panel: phase space of the Hamiltonian ${\overline{\mathcal{H}}}_0$ for the Dione–Helene system. Right panel: shift of the trailing libration center (dashed line) in our model with respect to the simple co-orbital motion (full line).

Fig. 2 Schematic orientation of the orbits for the Trojan satellites in the ACR condition, assuming λ 3 = 0° and ψ 2 = 0°.

Fig. 3 Topology of the inner evection resonance for different mass ratios. Panels a and b approximately correspond to the Dione–Helene pair. In all cases, a1 = a2 = 0.00326 au (488 000 km) and Δϖ = ϖ2 −ϖ1 = 60°.

Fig. 4 Topology of the outer evection resonance for m2 ≪ m1. Panels a andb: behavior of m1 at two different distances from the planet (in all cases a1 = a2). In panel a, two libration islands at ψ1 = 0° and 180° are connected by an unstable point that is shifted to positive values of e1 cos ψ 1. In panel b only a kinematic libration around ψ1 = 180° is observed. Panel c: behavior of m2 at the same distance from the planet as in panel a. In this case, the topology is the same, but the libration islands are rotated by −Δϖ. In all cases, Δϖ = 60°.

2.2 Application to the Saturn system

We applied our model to a physical system consisting of Saturn with m₀ = 2.858 × 10⁻⁴ M_⊙, R₀ = 60 268 km, J₂ = 1.6298 × 10⁻², J₄ = 0, a distance to the Sun a₃ = 9.537 au, and assuming e₃ = 0 (we recall that the eccentricity of Saturn is 0.054). We considered satellite masses on the order of those of Tethys and Dione or lower (≤ 6–10 × 10²⁰ kg), and varied the mass ratio m₂∕m₁ to describe the dynamics of the model. Hereafter, we refer to the angles ψ _i as the evectionangles.

Two different resonant regimes can be distinguished. The first occurs at very close distances from the central body and is driven by the perturbation of the planet oblateness. This regime corresponds to the inner evection resonance, and its topology displays twosymmetric libration islands that resemble a second-order Andoyer Hamiltonian. For low mass ratios m₂ ∕m₁, the libration center of the trailing Trojan, m₁, is located atψ₁ = ±90° (Fig. 3a). This topology is similar to the case of a single satellite (i.e., m₂ = 0). On the other hand, the behavior of the leading Trojan m₂ is exactly the same, but the libration centers are shifted to ± 90° −Δϖ (Fig. 3b). When m₂ = m₁, the topology is preserved, but the libration centers are shifted to ± 90° + Δϖ∕2 for the trailing and to to ± 90° −Δϖ∕2 for the leading Trojan (Figs. 3c and d). It is worth noting that the level curves in were computed for the same semi-major axes (i.e., the same values of Z₁), implying that the distance at which the inner evection resonance occurs is roughly independent of the mass ratio.

The second resonant regime occurs at large distances from the planet and is driven by the solar perturbation. This regime corresponds to the outer (or classical) evection resonance. The basic topology shows two asymmetric libration islands enclosed by a separatrix generated by a single saddle point close to the origin (Fig. 4a). The asymmetry is caused by the presence in the disturbing function of the odd-degree terms in eccentricity, as previously shown by Frouard et al. (2010), in particular the term e_i cosψ_i, which is associated with the variation harmonic in the classical lunar theory. For low m₂ ∕m₁ ratios, the libration center of the trailing Trojan is located at either ψ₁ = 180° or 0°, and this is again similar to the evection topology in the case of a single satellite. It is worth noting that for closer distances to the planet, the topology is not strictly resonant because there is no separatrix (Fig. 4b). In this case, we observe a kinematic libration around ψ ₁ = 180° forced by the odd-degree terms in e. On the other hand, since the unstable equilibrium that generates the separatrix is always shifted from the origin, a resonant orbit around ψ ₁ = 180° but very close to the separatrix may eventually display a circulating angle ψ₁. Figure 4c shows the topology of the leading Trojan, which is the same as in Fig. 4a, but rotated by −Δϖ.

We recall that in the case of a single satellite, the precession rate caused by the oblateness of Saturn is ${\stackrel{.}{\varpi }}_{\text{obl}}\approx 3n_1{J}_2{R}_0^2/(2a_1^2)$ (Roy 1978), while the precession rate caused by solar perturbation is ${\stackrel{.}{\varpi }}_{\text{sol}}\approx 3G{m}_3/(4n_1{a}_3^3)$ (Innanen et al. 1997). Then, according to Li & Christou (2016), the critical semi-major axis at which one precession rate dominates is $$a_{\text{crit}}={\left(2\frac{m_0}{m_3}{J}_2{R}_0^2{a}_3^3\right)}^{1/5}~0.0167\text{au},$$(23)

which defines the transition limit between the inner and outer evections in the Saturn system.

The locations of the libration centers of the inner and outer evection resonances, computed for the parameters of the Dione–Helene system (m₂∕m₁ ≃ 2 × 10⁻⁵), are shown in Fig. 5. These locations were obtained by searching for the conditions in the a, e plane that make $\partial \overline{\mathcal{H}}/\partial W_1=0$. The outer evection shows two branches, one corresponding to ψ₁ = 180°, which is always stable (blue line), and the other corresponding to ψ ₁ = 0°, which splits into the stable (blue line) and the unstable (red line) equilibria.

Figure 5 also shows the position of the prograde moons of Saturn. We note that all the moons are located inside of the location of the outer evection resonance, including the group of retrograde moons (not shown), whose outermost limit is indicated by the vertical dashed line. We also note that the regular moons, that is, those with e≲0.02 and i≲1°, are clustered inside of the location of the inner evection resonance (a≲0.003 au). In particular, all the moons involved in co-orbital configurations (Janus-Epimetheus, Tethys, and Dione) are in this group. The only exception among the regular moons is Rhea, which is located slightly outside of the inner evection, but has no co-orbitals.

In the case of equal masses, the topology of the outer evection resonance is preserved, but the libration centers of both the trailing and leading Trojans are rotated by ±Δϖ∕2, respectively (e.g., Fig. 6a). At variance with the inner evection resonance, the existence existence of the libration islands in this case is very sensitive to the value of Δϖ, and they quickly disappear if the Trojan pair moves away from the ACR solution toward increasing Δ ϖ (Fig. 6b). Another major difference with respect to the inner evection is that the distance from the planet at which the resonant islands appear (more precisely, where the unstable equilibrium appears) is strongly dependent on the mass ratio, as shown in Fig. 6c (full line). This figure also shows how the libration center moves away from ψ ₁ = 0 as the mass ratio tends to 1 (dashed line).

Fig. 5 Location of the inner and outer evections (full lines) for the masses and eccentricities of the Dione–Helene system. The stable branches (related to the libration centers) are shown in blue, and the unstable branch (related to the saddle point of the outer evection) is plotted in red. The full dots indicate the current position of the prograde satellites of Saturn, and the vertical dashed line indicates the outermost limit of the irregular satellites. The letters specify the locations of Dione (D), Tethys (T), and Rhea (R), respectively. For the outer evection, the stable branch at ψ 1 = 180° and e1 ≲ 0.14 is related to a kinematic libration (Fig. 4b).

3 Numerical study

In this section we present a numerical survey of the stability of the Saturn satellites, with particular focus on the occurrence of the evection resonance in both the cases of a single satellite and a Trojan pair. Our numerical model takes into account the perturbations from the Sun as well as from the oblateness of Saturn by including the J₂ zonal harmonic. The perturbations from other satellites in the Saturn system are not taken into account. The satellite orbits are described by the orbital elements a, e, λ, and ϖ. Hereafter, we maintain the subscripts 1,2 for the Trojan pair (trailing and leading, respectively), and we use the subscripts ⊙ and S for the Sun and Saturn, respectively.

The Newton equations of motion are numerically solved using a Burlisch–Stoer integrator with adaptive step size, which was modified to independently monitor the error in each variable; this imposes a relative precision better than 10⁻¹³. We stopped the integrations when the mutual distance between Saturn and the satellite was smaller than the sum of their radii. The total integration time was set to encompass several periods of the orbital secular variations.

The stability was evaluated through the use of some chaos and resonance indicators, and we explored different grids of initial conditions in the space of orbital elements. In particular, for each initial condition of the grid, we computed the value of the indicator for the mean exponential growth of nearby orbits (MEGNO), ${\langle Y\rangle }_i$, i = 1, 2 (Cincotta & Simó 2000; Maffione et al. 2011). This indicator is very efficient to quickly distinguish regular (${\langle Y\rangle }_i~2$) from chaotic orbits (${\langle Y\rangle }_i\gg 2$), but it does not provide a detailed representation of the structure of a resonance.

To better address the structure of the resonances, we also computed the amplitude of maximum variation of the orbital eccentricity of the satellites during the integrations, that is, $\delta e_i=\left(e_{i,\text{max}}-e_{i,\text{min}}\right)/2$, i = 1, 2. This indicator has proven to be an extremely useful tool for mapping the resonant structure in N-body problems (e.g., Ramos et al. 2015). We recall that δe_i is not a direct measure of chaotic motion, but abrupt changes in δe_i are often tracers for the presence of resonances. Therefore, regions with large variations in δ e_i are more sensitive to perturbations and are very likely chaotic. Frouard et al. (2010) used a similar indicator in their study of the evection resonance in Jupiter satellites, defined as $(e_{\text{max}}^{{\tau }_a}-e_{\text{max}}^{{\tau }_b})/e_{\text{max}}^{{\tau }_a}$, where the total integration time, τ, was divided into two consecutive samples τ_a and τ_b. Nevertheless, we realized that this indicator is quite sensitive to the value of τ and may artificially create some structures that are not related with actual resonances.

A third indicator that we used is the amplitude of oscillation of the evection angle, defined as $\delta {\psi }_i=\left|{\psi }_{i,\text{max}}-{\psi }_{i,\text{min}}\right|$. This indicator has a minimum at the center of the resonance, or more generally, at the condition ${\stackrel{.}{\psi }}_i={\stackrel{.}{\lambda }}_{\odot }-{\stackrel{.}{\varpi }}_i~0$.

Fig. 6 Behavior of the outer evection resonance for equal masses of the Trojan pair. Panel a: topology of m1 at a distance of a1 = a2 = 0.39 au (5.88 × 107 km), and assuming a periastron difference Δϖ = 60° (i.e., atthe ACR solution). Panel b: same as panel a, but for a periastron difference Δ ϖ = 70°. Panel c: distance from the planet at which the evection resonance appears (full line), and the angle of the resonance center (dashed line) as a function of the mass ratio.

Fig. 7 Dynamical maps in the (a, e)-plane for the regions of the inner evection (top row) and the outer evection (bottom row). The MEGNO, δ e, and δ ψ indicators are displayed in the left, middle, and right columns, respectively. Angles are in degrees. White areas correspond to unstable orbits. The gray areas in the δ ψ maps indicate non-librating orbits. The grids correspond to osculating orbital elements with initial ψ = 90° for the inner evection, and initial ψ = 0° for the outer evection. To allow comparison with previous works, the upper horizontal axes give the distance to the planet in convenient units: planetary radii (RS) for the inner evection, and Hill radii (RHill) for the outer evection.

3.1 Evection for a single moon

Figure 7 shows a set of dynamical maps in the (a, e)-plane for the inner (top row, ψ = 90°) and the outer (bottom row, ψ = 0°) evections. The maps display the values of the three indicators described in Sect. 3: $\langle Y\rangle$ (left column), δ e (middle column), and δψ (right column). We drop here the subindex i for simplicity.

The inner evection region is quite regular, and both δe and δψ allow identifying the exact location of the resonance. For eccentricities lower than ~ 0.6, this location coincides with the eccentricity predicted by the analytical model (cf. Fig. 5). At higher eccentricities, the resonance location moves toward slightly larger semi-major axes, a behavior that the analytical model does not reproduce because the ${\overline{\mathcal{J}}}_i$ series is truncated early. The gray areas in the map of δψ correspond to conditions where the angle ψ circulates. White areas correspond to orbits that either collide with Saturn (a ≤ R_S) or escape from the system (a > a_S) during the integration time-span.

The outer evection resonance, on the other hand, occurs inside a stable island surrounded by a thin layer of chaotic motion (red in the map of $\langle Y\rangle$). The map of δ ψ shows that the libration amplitude around ψ = 0° is never larger than ~60°, in good agreement with the analytical model (cf. Fig. 4a). The map of δ e also shows a stable narrow region at low eccentricities (blue) that is related to the kinematic librations induced by the forced eccentricity mode. Nevertheless, we note that the resonant island occurs closer to the planet (a ~ 0.19 au) than predicted by the analytical model (a ~ 0.26 au). This is due in part to the fact that the grids in Fig. 7 are given in terms of the osculating initial semi-major axis, while in terms of the mean semi-major axis, the island shifts to higher values, as shown in Fig. 8. Another reason might be the early truncation of the disturbing function Eq. (10) and that the eccentricity of the Sun is not considered in the model. A similar behavior has been identified by Frouard et al. (2010).

Figure 8 shows the dynamical maps of δ ψ in the (a, ψ)-plane for the outer evection. In this case, we have initially e = 0.4 and λ_⊙ = 0°, which means that varying the initial ψ is equivalent to varying the initial ϖ. We note that the libration regions around 0° or 180° occur at different osculating semi-major axes (Fig. 8, left), but this dependence disappears when the mean semimajor axis is considered (Fig. 8, right).

Since the inner evection is driven by the planet oblateness, it is expected that the planet obliquity plays some role in the behavior of the resonance. In Fig. 9 we address the effect of the obliquity of Saturn (ϵ _S = 26.73°) and its orbital eccentricity (e_S = 0.054) on the inner evection. We note that a non-zero obliquity shifts the location of the resonance to slightly smaller distances (from ~ 8 to ~ 7 R_S), although the structure of the equilibrium family is preserved. On the other hand, the eccentricity of the Saturn orbit (or equivalently, the Sun) does not produce any significant changes.

Fig. 8 Dynamical maps of the δψ indicator in the (a, ψ)-plane for the outer evection, with initial e = 0.4. The leftpanel corresponds to a grid of osculating elements, while the right panel is the same grid in terms of mean elements. Blue regions represent the conditions for which δ ψ ≲ 30°. Gray regionscorrespond to larger libration amplitudes or circulation of ψ. White regions are unstable orbits. The resonance occurs at a mean semi-major axis ā = 0.205 au. The upper horizontal axis is the distance in terms of RHill.

Fig. 9 Value of δe as a function of the distance from the planet in the region of the inner evection and for e = 0. The curve in violet corresponds to the case of an obliquity ϵS = 0° and orbital eccentricity eS = 0° for Saturn. The green dots correspond to the case of ϵS = 26.73°, either with eS = 0° or eS = 0.054°. The minimum inside the peak of the curves marks the center of the evection resonance.

Fig. 10 Dynamical map of the δψ1 indicator in the osculating (a, e)-plane for the inner evection, setting ψ1 = 120° initially. Case of Trojan satellites with equal masses. The resonance occurs at the same location as in the case of a single satellite (compare with Fig. 7). The upper horizontal axis gives the distance in RS. Angles are in degrees.

3.2 Evection for Trojan moons

In this section, we study the evolution of two co-orbital moons in Trojan configuration initially near the equilibrium of the ACR condition, that is, σ ≡Δλ = 60°, Δ ϖ= 60°, e₂ = e₁ (Giuppone et al. 2010; Robutel & Pousse 2013). We consider two cases: m₂ ∕m₂ = 1 and m₂ ∕m₁ = 10⁻⁴. All the simulations initially have λ _⊙ = 0°. Figure 2 shows a schematic orbital diagram when ψ₂ = 0°.

Figure 10 shows the dynamical map of the δ ψ _i indicator for a system of two Trojan moons with equal masses in the region of the inner evection. The indicator corresponds to those of the trailing Trojan (i = 1), but the situation is identical for the leading Trojan. A comparison with Fig. 7 shows that the evection resonance occurs at almost the same distance from the planet as for a single moon (a ~ 0.003 au at e~ 0), implying that the location of the inner evection is independent of the mass ratio of the Trojans. This agrees well with the predictions of the analytical model (Fig. 3).

On the other hand, the outer evection in the case of Trojan moons cannot be studied because the Trojan configuration does not survive at large distances from Saturn. This is shown in Fig. 11 for moons of equal masses. The white regions correspond to the conditions where an initial Trojan ACR configuration is destroyed, and the co-orbital pair does not survive as such. Stable Trojan configurations exist up to distances of ~ 0.06 au only for almost circular orbits. At distances ≲0.02 au (blue regions), the Trojan ACR configuration is stable even for high eccentricities, except at the locations of the evection resonance and its higher harmonics (${\stackrel{.}{\varpi }}_i~{\stackrel{.}{\lambda }}_{\odot }$; ${\stackrel{.}{\varpi }}_i~2{\stackrel{.}{\lambda }}_{\odot };2{\stackrel{.}{\lambda }}_{\odot }$; ${\stackrel{.}{\varpi }}_i~3{\stackrel{.}{\lambda }}_{\odot }$). The same qualitative results are obtained for a Trojan pair with In the region where the outer evection resonance occurs for single moons (a ~ 0.2 au; see Fig. 7), Trojan configurations are not sustainable for more than 100 yr. As a consequence, a pair of moons initially inside the evection island and in a Trojan ACR configuration quickly ends up as a non-Trojan pair with only one of the components trapped in the evection resonance. This is shown in Fig. 12 for a pair of equal-mass moons. The top panel shows the case where ψ ₁ librates while ψ₂ circulates, and the bottom panel shows the opposite case. Neither of these cases show σ librating around ± 60°, so Trojan motion no longer exists. Simultaneous libration of ψ₁, ψ₂ is never found in the outer evection either. Figure 13 shows the same evolutions as Fig. 12, but in the space k_i = e_i cosψ_i, h_i = e_i sinψ_i.

The libration of either ψ₁ or ψ ₂ is related to the initial longitudes of the periastrons relative to the longitude of the Sun (e.g., Fig. 2). The behavior of the former Trojan companion that becomes capture in the evection resonance is very similar to that of a single moon. In the (a, e)-plane, the resonance occurs inside an island surrounded by a layer of chaotic motion, as shown in Fig. 14. We note, however, that the island corresponding to the libration of ψ ₁ (left panel) is smaller than that corresponding to the libration of ψ₂ (right panel), indicating that even in the case of equal masses, the configuration may not be symmetric. Neither configuration is long-term stable, and at some point, the moon where ψ _i circulates is ejected or hits the other moon or the planet. Nevertheless, configurations that lead to the libration of ψ ₁ appear to be more stable (up to 10⁶ yr) than those leading to the libration of ψ₂ (up to ~10³ yr; see Fig. 13).

Fig. 11 Dynamical maps in the osculating (a, e)-plane displaying the amplitude δe1 of the eccentricity variation (top panel), the amplitude of the σ libration (middle panel), and the amplitude of the Δϖ libration (bottom panel) for a system of two equal-mass moons initially in a Trojan ACR configuration. Angles are in degrees. Trojan motion is unstable in the white regions in less than 5 × 104 yr, and it is most stable in the blue regions. The stripes of less stable Trojan motion observed for a < 0.03 au correspond to the locations of the evection resonance and its higher harmonics.

4 Tidal evolution

The dynamical evolution of the Saturn system under the effects of tidal forces is relevant for constrainingthe past evolution of the moons. In particular, the proximity of the inner evection to the planet and to the few co-orbital satellites of the system raises the question about the possible effect of this resonance on a pair of Trojan moons that evolve by tidal effects.

Short and long tidal migration rates of the Saturn moons have been explored in a wide variety of scenarios by Ćuk et al. (2016). The authors concluded that either the moons are significantly younger than the planet (formed only about 100 Myr ago), or their tidal evolution must be extremely slow (implying a dissipation factor for Saturn of Q > 80 000).

They also concluded that the inner evection resonance is located in a chaotic region, possibly because several harmonics of their own resonance overlap, although we showed here that this is not the case for isolated moons or Trojan pairs (see Figs. 7 and 10). Ćuk et al. (2016) proposed a mechanism to form a disk of proto-moons that could produce a favorable environment for moons to be created, which increases the possibility that some of them might harbor Trojan companions.

The dissipation factor of Saturn has traditionally been estimated to be Q ~ 18 000, which gives a very low tidal migration rate. With this value, Enceladus, Tethys, Dione, and Rhea are not expected to have migrated significantly over the age of the solar system; they were probably formed at their current locations. However, based on an analysis of astrometric positions from old photographic plates and images of the Hubble Space Telescope, Lainey et al. (2012) revealed that the Saturn moons did migrate over the last century. They estimated a dissipation factor ten times smaller, Q ~ 1 682 ± 540. Based on this result, Crida & Charnoz (2014) suggested that the inner moons, up to Rhea, could have been formed from the primordial spreading ring of Saturn, and in this scenario, it may be expected that many primordial Trojan pairs have crossed the evection resonance.

Here, we study the effect of tidal forces in the evolution of hypothetical Trojan moons by incorporating a weak-friction equilibrium tidal model into our N-body code. The tidal forces were modeled using the closed formulas of Mignard (1979) and following the approach of Rodríguez et al. (2013) for co-orbital companions. For Saturn, we assumed a Love number k₂ = 0.39 and Q = 1 682 (Lainey et al. 2012), and for the satellites, we assumed k₂ = 0.05 and 20 ≤ Q ≤ 200 (Ćuk et al. 2016).

Numerical integrations using the nominal values of the system parameters are impracticable. Therefore, we scaled the total time-span and the tidal parameters down by the same factor α > 1 to reduce the total integration time, as suggested by Zoppetti et al. (2018). We considered values of α low enough to keep the accelerated tidal evolution adiabatic with respect to the resonant and secular gravitational timescales. The results shown here correspond to α = 200, but we performed some other numerical test using lower values to check the reliability of the results.

We know that in the adiabatic regime, the system evolves along the family of periodic orbits, thus we ran some simulations with a single satellite. We considered a fictitious Dione-like moon that was captured in the evection resonance at a = 0.00326 au, with ψ librating around ± 90°, as predicted by the analytical model. When we increased the value of α, the libration amplitude increased and the satellite was able to cross the evection resonance without any effect on it.

Then we ran experiments considering a pair of Trojan moons with mass ratios of m₂ ∕m₁ = 2 × 10⁻⁵, m₂ ∕m₁ = 1, and m₂ ∕m₁ = 2 × 10⁵. We recall that the first ratio approximately corresponds to the Dione–Helene system. In all cases, the initial semi-major axis was set to a = 0.00321 au, slightly inside of the evection resonance, and the initial eccentricity was set to the current values of Dione and Helene (0.0022). Typical results are shown in Fig. 15, where the evolution of m₁ is identify in magenta and that of m₂ in green. The resonant capture is observed in the bottom panels, where both evection angles ψ₁, ψ₂ start to librate. For very low or very high mass ratios, the most massive component is always captured around 90° or 270°, while the less massive componentis captured around ± 60° of those values. For equal masses, the components are captured around 60° (or 140°) and 120° (or 300°), respectively,as predicted by the analytical model. The evolution of the Trojan pair inside the evection resonance implies that the eccentricities start to increase, moving along the family of periodic orbits with oscillations of higher amplitude for the less massive component. We also note that the ACR equilibrium configuration is preserved during the evolution inside the evection resonance. For low mass ratios of the Trojan pair, excitation of the eccentricities contributes to destabilize the system and may lead to the loss of the less massive component, either by hitting the planet or the other component.

The results are similar over the whole range of values of Q we tested for the satellites. On the other hand, when we considered lower values of Q for Saturn (Q < 1 000), which implies abandoning the adiabatic regime, the evection resonance causes a small jump in the eccentricities of the Trojan satellites that breaks the libration of Δϖ, causing the system to depart from the ACR solution. After this, the less massive satellite is captured with with ψ _i librating around ± 90°, the angle Δϖ starts to circulate, and σ → 180°, leading the pair into a horseshoe regime. The system survives under this condition until a collision occurs after ~ 10⁶ yr or sooner.

Fig. 12 Evolution of the resonant angles σ (cyan lines) and ψi (red/blue lines) for a pair of moons with equal masses, initially placed in a Trojan ACR configuration (σ = 60°, Δ ϖ = 60°) in the region of the outer evection (a1 = a2 ~ 0.2 au, e1 = e2= 0.4). Top panel: case ψ1 = 0° and a1 = a2 = 0.1825 au initially. Bottom panel: case ψ2 = 0° and a1 = a2 = 0.1893 au initially.

Fig. 13 Same as Fig. 12, but in the ki, hi space (see text) and for the full time span of the simulation. The top panels correspond to the top panel of Fig. 12; the system survives for at least 106 yr. The bottom panels correspond to the bottom panel of Fig. 12; the system survives for ~ 103 yr. In both cases, λ⊙ = 0° initially.

Fig. 14 Dynamical maps of δψi in the osculating (a, e)-plane for a Trojan pair of equal masses, initially placed in the ACR configuration. Initial conditions for the angles are the same as in Fig. 13. Blue islands correspond to libration of ψ 1 (left panel) or ψ2 (right panel). Gray regions correspond to circulation of ψi, and white regions correspond to unstable orbits (compare to Fig. 7). The upper horizontal axis is the distance in terms of RHill.

Fig. 15 Adiabatic evolution of pair of Trojan satellites that migrate as a result of tidal effects, using Q = 100. Each column of panels corresponds to a different mass ratio m2∕m1: 2 × 10−5 (left panels),1 (middle panels), and 2 × 105 (right panels). Eccentricity is shown in the top row, σ and Δ ϖ in the middle row, and ψ1, ψ2 in the bottom row. Magenta and green identify identify the trailing and leading Trojan, respectively, except in the middle panels, where they identify σ and Δ ϖ,

5 Conclusions

We have addressed the dynamics of the evection resonance in the case of co-orbital motion, and focused on Trojan pairs in the ACR configuration (Δλ = 60°, Δϖ = 60°). We applied an expansion of the averaged Hamiltonian describing the motion of two bodies in Trojan orbits, perturbed by a third distant body. The expansion combines the model developed by Frouard et al. (2010) for the evection of single bodies and the model described by Robutel & Pousse (2013) for the motion of Trojan bodies. This model can be applied for Trojan satellites in the solar system as well as for circumstellar Trojan planets in binary systems. Here we consider the Saturn system. We also performed a series of numerical simulations, both to validate the predictions of the model and to better assess the actual dynamics of Trojan satellites around Saturn. Our conclusions are summarized below.

• Our model addressed the two types of evection resonance that may appear: the inner evection resonance, close to the planet and related to the precession of the periastron forced by the oblateness of the central body, and the outer evection resonance, far away from the planet and related to the precession of the periastron forced by the distant perturber.

• The global structure of the phase space of the evection resonances for Trojan satellites is similar to that of a single satellite, differing in that the libration centers are displaced from their standard positions (± 90° for the inner and 0°, 180° for the outer evection) by an angle that depends on the periastron difference ϖ₂ −ϖ₁ and on the mass ratio m₂∕m₁ of the Trojan pair.

• The distance from the planet at which the inner evection resonance occurs is independent of the mass ratio of the Trojanpair. The resonance only moves toward slightly smaller distances when the obliquity of the planet is taken into account.

• In the Saturn system, the evolution of Trojan moons into the inner evection resonance, located at ~8 R_S, is possible. Capture in this resonance can be driven by adiabatic tidal migration, exciting the eccentricities of the two Trojan components. This in turn may destabilize the resonant system, leading to the loss of one of the components.

• The location of the outer evection resonance depends, in principle, on the mass ratio of the Trojan pair, moving to longer distances as m₂∕m₁ → 1. For the Saturn system, this would occur at distances > 0.4 R_Hill.

• Nevertheless, the outer evection in the Saturn system cannot exist at all for Trojan moons because Trojan configurations are strongly unstable at distances from Saturn longer than ~0.15 R_Hill. A pair of moons forced into a Trojan configuration in the region of the outer evection always ends in a non-co-orbital configuration, with one of the moons captured in a stable orbit inside the evection resonance, while the other moon becomes unstable and is eventually lost.

In the Saturn system, the moons that are closer to the current location of the inner evection are Tethys, Dione, and Rhea. Of these, only the first two have Trojan companions and are located at smaller distances from Saturn than the evection. Rhea has no Trojans andis located slightly beyond the evection. If Rhea had Trojan companions in the past, it may have lost them when it crossed the evection resonance. Thus, it would be interesting to study the possibility that the craters on the surface of Rhea were produced by collisions with former Trojan companions that were destabilized by the evection. This study might be extended even for other regular satellites such as Titan, Hyperion, and Iapetus.

The high eccentricity of Titan, Hyperion, and Iapetus might be related to the excitation caused by the crossing of the inner evection resonance. The eccentricity of Rhea, on the other hand, may have been excited by the evection and may have later been damped to the current value by tidal dissipation.

On the other hand, the inner evection resonance moved in the past from its current location as a result of changes in the shape or obliquity of Saturn and because Saturn was closer to the Sun than it is today (e.g., Nesvorný & Morbidelli 2012). Then, Tethys and Dione might have crossed the resonance inthe past and preserved their Trojan companions. This might impose constraints on the dynamical processes that led to this type of evolution.

As for the outer evection, none of the current Saturn moons (not even the irregular moons) reach the current domain of this resonance. However, the resonance may have been relevant for the early evolution of the moons, when Saturn was closer to the Sun.

Acknowledgements

The authors acknowledge useful discussions with C. Beaugé that helped us to improve the speed of the numerical integrations with the Bulirsch–Stoer method. We also wish to thank the anonymous referee for their comments. This work has been supported by the National Council of Research of Argentina (CONICET) and the Rio de Janeiro State Science Foundation (FAPERJ), in the framework of a FAPERJ/CONICET bilateral project, and by the Brazilian National Council of Research (CNPq). Computations were performed at the BlaFis cluster of the University of Aveiro, Portugal, and at the SDumont cluster of the Brazilian System of HPC (SINAPAD).

Appendix A Coefficients of $\overline{\mathcal{H}}$

We provide here the explicit expressions of the coefficient functions of the Hamiltonian Eq. (20) in terms of non-canonical orbital elements. $$A_{0i}=-G{m}_3{m}_i{a}_i^2\left(1-{\nu }_i\right)\frac{1}{4},$$ $$A_{1i}=G{m}_3{m}_i{a}_i^3\frac{5}{16},$$ $$C_0=G{m}_3\overline{\nu }{a}_1{a}_2\frac{1}{2}\cos 2{\theta }_1,$$ $$C_1=G{m}_3\overline{\nu }{a}_1{a}_2\frac{1}{8}\cos 4{\theta }_1,$$ $$C_2=G{m}_3\overline{\nu }{a}_1{a}_2\frac{9}{8},$$ $$C_3=G{m}_3\overline{\nu }{a}_1{a}_2\frac{1}{8}\sin 4{\theta }_1,$$ $$C_4=G{m}_3\overline{\nu }{a}_1{a}_2\frac{3}{16}\sin 2{\theta }_1,$$ $$D_i=-G{m}_0{m}_i\frac{1}{2}{J}_2\frac{R_0^2}{a_i^3},$$ $${\text{Δ}}_{12}={\left(a_1^2+a_2^2-2a_1{a}_2\cos 2{\theta }_1\right)}^{1/2},$$ $$B_0=-{\sum }_{i=1}^2\frac{{\beta }_i{\mu }_i}{2a_i}\frac{2{\overline{\nu }}^2}{{\mu }_1{\mu }_2}-\frac{G{m}_1{m}_2}{{\text{Δ}}_{12}}+\overline{\nu }\left(1+\frac{2\overline{\nu }}{m_0}\right)\sqrt{\frac{{\beta }_1{\beta }_2}{a_1{a}_2}}\cos 2{\theta }_1,$$ $$\begin{array}{ll}{B}_1=\hfill & G{m}_1{m}_2\left[\frac{a_1{a}_2(a_1^2+a_2^2)\cos 2{\theta }_1}{2{\text{Δ}}_{12}^5}-\frac{a_1^2{a}_2^2\left(13-5\cos 4{\theta }_1\right)}{8{\text{Δ}}_{12}^5}\right]\hfill \\ \hfill & -\overline{\nu }\left(1+\frac{2\overline{\nu }}{m_0}\right)\sqrt{\frac{{\beta }_1{\beta }_2}{a_1{a}_2}}\frac{1}{2}\cos 2{\theta }_1,\hfill \end{array}$$ $$\begin{array}{ll}{B}_2=\hfill & -\text{}G{m}_1{m}_2[\frac{a_1{a}_2(a_1^2+a_2^2)\cos 4{\theta }_1}{{\text{Δ}}_{12}^5}\hfill \\ \hfill & +\frac{a_1^2{a}_2^2\left(\cos 6{\theta }_1-17\cos 2{\theta }_1\right)}{8{\text{Δ}}_{12}^5}]\hfill \\ \hfill & +\overline{\nu }\left(1+\frac{2\overline{\nu }}{m_0}\right)\sqrt{\frac{{\beta }_1{\beta }_2}{a_1{a}_2}}\cos 4{\theta }_1,\hfill \end{array}$$ $$\begin{array}{ll}{B}_3=\hfill & -\text{}G{m}_1{m}_2[\frac{a_1{a}_2(a_1^2+a_2^2)\sin 4{\theta }_1}{{\text{Δ}}_{12}^5}\hfill \\ \hfill & +\frac{a_1^2{a}_2^2\left(\sin 6{\theta }_1-35\sin 2{\theta }_1\right)}{8{\text{Δ}}_{12}^5}]\hfill \\ \hfill & +\overline{\nu }\left(1+\frac{2\overline{\nu }}{m_0}\right)\sqrt{\frac{{\beta }_1{\beta }_2}{a_1{a}_2}}\sin 4{\theta }_1,\hfill \end{array}$$

References

All Figures

	Fig. 1 Left panel: phase space of the Hamiltonian ${\overline{\mathcal{H}}}_0$ for the Dione–Helene system. Right panel: shift of the trailing libration center (dashed line) in our model with respect to the simple co-orbital motion (full line).
In the text

	Fig. 2 Schematic orientation of the orbits for the Trojan satellites in the ACR condition, assuming λ 3 = 0° and ψ 2 = 0°.
In the text

	Fig. 3 Topology of the inner evection resonance for different mass ratios. Panels a and b approximately correspond to the Dione–Helene pair. In all cases, a1 = a2 = 0.00326 au (488 000 km) and Δϖ = ϖ2 −ϖ1 = 60°.
In the text

Fig. 4 Topology of the outer evection resonance for m2 ≪ m1. Panels a andb: behavior of m1 at two different distances from the planet (in all cases a1 = a2). In panel a, two libration islands at ψ1 = 0° and 180° are connected by an unstable point that is shifted to positive values of e1 cos ψ 1. In panel b only a kinematic libration around ψ1 = 180° is observed. Panel c: behavior of m2 at the same distance from the planet as in panel a. In this case, the topology is the same, but the libration islands are rotated by −Δϖ. In all cases, Δϖ = 60°.

In the text

Fig. 5 Location of the inner and outer evections (full lines) for the masses and eccentricities of the Dione–Helene system. The stable branches (related to the libration centers) are shown in blue, and the unstable branch (related to the saddle point of the outer evection) is plotted in red. The full dots indicate the current position of the prograde satellites of Saturn, and the vertical dashed line indicates the outermost limit of the irregular satellites. The letters specify the locations of Dione (D), Tethys (T), and Rhea (R), respectively. For the outer evection, the stable branch at ψ 1 = 180° and e1 ≲ 0.14 is related to a kinematic libration (Fig. 4b).

In the text

Fig. 6 Behavior of the outer evection resonance for equal masses of the Trojan pair. Panel a: topology of m1 at a distance of a1 = a2 = 0.39 au (5.88 × 107 km), and assuming a periastron difference Δϖ = 60° (i.e., atthe ACR solution). Panel b: same as panel a, but for a periastron difference Δ ϖ = 70°. Panel c: distance from the planet at which the evection resonance appears (full line), and the angle of the resonance center (dashed line) as a function of the mass ratio.

In the text

Fig. 7 Dynamical maps in the (a, e)-plane for the regions of the inner evection (top row) and the outer evection (bottom row). The MEGNO, δ e, and δ ψ indicators are displayed in the left, middle, and right columns, respectively. Angles are in degrees. White areas correspond to unstable orbits. The gray areas in the δ ψ maps indicate non-librating orbits. The grids correspond to osculating orbital elements with initial ψ = 90° for the inner evection, and initial ψ = 0° for the outer evection. To allow comparison with previous works, the upper horizontal axes give the distance to the planet in convenient units: planetary radii (RS) for the inner evection, and Hill radii (RHill) for the outer evection.

In the text

Fig. 8 Dynamical maps of the δψ indicator in the (a, ψ)-plane for the outer evection, with initial e = 0.4. The leftpanel corresponds to a grid of osculating elements, while the right panel is the same grid in terms of mean elements. Blue regions represent the conditions for which δ ψ ≲ 30°. Gray regionscorrespond to larger libration amplitudes or circulation of ψ. White regions are unstable orbits. The resonance occurs at a mean semi-major axis ā = 0.205 au. The upper horizontal axis is the distance in terms of RHill.

In the text

Fig. 9 Value of δe as a function of the distance from the planet in the region of the inner evection and for e = 0. The curve in violet corresponds to the case of an obliquity ϵS = 0° and orbital eccentricity eS = 0° for Saturn. The green dots correspond to the case of ϵS = 26.73°, either with eS = 0° or eS = 0.054°. The minimum inside the peak of the curves marks the center of the evection resonance.

In the text

	Fig. 10 Dynamical map of the δψ1 indicator in the osculating (a, e)-plane for the inner evection, setting ψ1 = 120° initially. Case of Trojan satellites with equal masses. The resonance occurs at the same location as in the case of a single satellite (compare with Fig. 7). The upper horizontal axis gives the distance in RS. Angles are in degrees.
In the text

Fig. 11 Dynamical maps in the osculating (a, e)-plane displaying the amplitude δe1 of the eccentricity variation (top panel), the amplitude of the σ libration (middle panel), and the amplitude of the Δϖ libration (bottom panel) for a system of two equal-mass moons initially in a Trojan ACR configuration. Angles are in degrees. Trojan motion is unstable in the white regions in less than 5 × 104 yr, and it is most stable in the blue regions. The stripes of less stable Trojan motion observed for a < 0.03 au correspond to the locations of the evection resonance and its higher harmonics.

In the text

	Fig. 12 Evolution of the resonant angles σ (cyan lines) and ψi (red/blue lines) for a pair of moons with equal masses, initially placed in a Trojan ACR configuration (σ = 60°, Δ ϖ = 60°) in the region of the outer evection (a1 = a2 ~ 0.2 au, e1 = e2= 0.4). Top panel: case ψ1 = 0° and a1 = a2 = 0.1825 au initially. Bottom panel: case ψ2 = 0° and a1 = a2 = 0.1893 au initially.
In the text

	Fig. 13 Same as Fig. 12, but in the ki, hi space (see text) and for the full time span of the simulation. The top panels correspond to the top panel of Fig. 12; the system survives for at least 106 yr. The bottom panels correspond to the bottom panel of Fig. 12; the system survives for ~ 103 yr. In both cases, λ⊙ = 0° initially.
In the text

Fig. 14 Dynamical maps of δψi in the osculating (a, e)-plane for a Trojan pair of equal masses, initially placed in the ACR configuration. Initial conditions for the angles are the same as in Fig. 13. Blue islands correspond to libration of ψ 1 (left panel) or ψ2 (right panel). Gray regions correspond to circulation of ψi, and white regions correspond to unstable orbits (compare to Fig. 7). The upper horizontal axis is the distance in terms of RHill.

In the text

Fig. 15 Adiabatic evolution of pair of Trojan satellites that migrate as a result of tidal effects, using Q = 100. Each column of panels corresponds to a different mass ratio m2∕m1: 2 × 10−5 (left panels),1 (middle panels), and 2 × 105 (right panels). Eccentricity is shown in the top row, σ and Δ ϖ in the middle row, and ψ1, ψ2 in the bottom row. Magenta and green identify identify the trailing and leading Trojan, respectively, except in the middle panels, where they identify σ and Δ ϖ,

In the text