11institutetext:  CIDMA, Departamento de Matemática, Universidade de Aveiro,
Campus de Santiago, 3810-193 Aveiro, Portugal 22institutetext: CFisUC,
Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal
33institutetext: ASD, IMCCE-CNRS UMR8028, Observatoire de Paris, PSL
Université, 77 Av. Denfert-Rochereau, 75014 Paris, France

# Eviction-like resonances for satellite orbits

Application to Phobos, the main satellite of Mars

Timothée Vaillant 1122 Alexandre C. M. Correia 2233

The motion of a satellite can experience secular resonances between the
precession frequencies of its orbit and the mean motion of the host planet
around the star. Some of these resonances can significantly modify the
eccentricity (evection resonance) and the inclination (eviction resonance) of
the satellite. In this paper, we study in detail the secular resonances that
can disturb the orbit of a satellite, in particular the eviction-like ones.
Although the inclination is always disturbed while crossing one eviction-like
resonance, capture can only occur when the semi-major axis is decreasing. This
is, for instance, the case of Phobos, the largest satellite of Mars, that will
cross some of these resonances in the future because its orbit is shrinking
owing to tidal effects. We estimate the impact of resonance crossing in the
orbit of the satellite, including the capture probabilities, as a function of
several parameters, such as the eccentricity and the inclination of the
satellite, and the obliquity of the planet. Finally, we use the method of the
frequency map analysis to study the resonant dynamics based on stability maps,
and we show that some of the secular resonances may overlap, which leads to
chaotic motion for the inclination of the satellite.

###### Key Words.:

celestial mechanics – planets and satellites: individual: Phobos – planets and
satellites: dynamical evolution and stability

##  1 Introduction

The orbit of a satellite can give some hints as to its origin. If the
satellite forms at the same time as its planet within an accretion disk, it is
expected that the satellite is located in the equatorial plane of the planet
in a nearly circular orbit (e.g., Batygin & Morbidelli 2020; Inderbitzi et al.
2020). If the satellite forms from a giant impact, the initial orbit is also
expected to be equatorial, but it can be very eccentric (e.g., Canup & Asphaug
2001; Canup 2005). Finally, if the satellite is a captured object, its orbit
is expected to be very eccentric as well and to present some inclination with
respect to the planet’s equator (e.g., Agnor & Hamilton 2006; Nesvorný et al.
2007).

Tidal dissipation inside the planet and the satellite can significantly modify
some orbital parameters, in particular the semi-major axis and the
eccentricity. Tidal dissipation can also lead to variations in the inclination
of the satellite, although these variations are weak and often negligible
(e.g., Lambeck 1979; Szeto 1983). Therefore, the inclination is not modified
by tides much and, contrary to the semi-major axis and the eccentricity, it
can be used to distinguish between the capture scenario and other formation
scenarios. Nevertheless, secular resonances can also modify the orbital
parameters of the satellite. For instance, in the case of Phobos, the main
satellite of Mars, Yoder (1982) noted that during its evolution, the satellite
crossed several resonances which have successively excited its eccentricity.
Therefore, the current value cannot be considered simply as the result of the
tidal dissipation of the primordial eccentricity. The evection resonance,
which corresponds to an interaction between the secular precession of the
pericenter of the satellite and the orbital mean motion of the planet, can
also lead to a variation in the eccentricity of the satellite (e.g., Touma &
Wisdom 1998). This effect is stronger when capture in the evection resonance
occurs, although this is only possible if the semi-major axis is increasing,
which is, for instance, the case for the early Moon.

There are also resonances that are able to modify the inclination of a
satellite, corresponding to an interaction between the precession of the node
of its orbit and the orbital mean motion of the planet. Touma & Wisdom (1998)
noted that the early Moon could also have been temporarily captured in a
resonance of this kind, which they called ”eviction” (just after it escapes
from the evection resonance). This temporary capture can lead to an increase
of around ten degrees in the inclination of the early Moon with respect to the
Earth’s equator, which may explain the current mutual inclination of the lunar
orbit of about five degrees. The capture of the Moon in the eviction resonance
is only possible if the semi-major axis is decreasing and if the eccentricity
is sufficiently high. Touma & Wisdom (1998) also noted that if the Moon
encounters the eviction resonance while the semi-major axis increases, it
would not be captured, but it could nevertheless lead to a one-time increase
in the inclination, provided that the eccentricity is high enough.

As the semi-major axis of Phobos is currently decreasing due to tides (e.g.,
Jacobson & Lainey 2014), it can also be captured in eviction-like resonances.
Indeed, Yokoyama (2002) showed that Phobos will encounter two of these
resonances in the near future, which may significantly increase its equatorial
inclination. Contrary to the original eviction resonance described by Touma &
Wisdom (1998), capture in these two resonances is possible for zero
eccentricity. In order to study these resonances, Yokoyama (2002) considered a
system of three bodies with the Sun, Mars, and a massless satellite, where the
orbit of Mars is circular. For this system, Yokoyama (2002) numerically
integrated the equations of the satellite obtained from an averaged
Hamiltonian, and noted that the capture probabilities in these resonances
strongly depend on the Martian obliquity. When capture does not occur, the
resonance crossing can still lead to a sudden variation in the inclination.
The capture probability in one resonance is almost certain if the obliquity of
Mars is higher than $20$ degrees. Similar results for this resonance are
observed in the case where the precession of the Martian equator, the
obliquity variations (Yokoyama 2002), or the planetary perturbations on Mars
are considered (Yokoyama et al. 2005). For the other resonance, which occurs
almost simultaneously with the evection resonance, Yokoyama (2002) numerically
observed that capture is possible with the resonant Hamiltonian, but becomes
impossible with the total Hamiltonian. Yokoyama (2002) then concluded that the
impossibility of the capture is due to the interaction between the two
resonances. Yokoyama et al. (2005) considered a more complete model including
the Mars’ eccentricity, and observed that it increases the interaction with
the evection resonance and its effects.

In this paper, we revisit the secular resonances between the orbit of the
satellite and the mean motion of the planet. We focus on those that are
suitable to excite the inclination of the satellite (eviction-like), since the
inclination can be used to put constrains on the formation scenarios. The aim
of the paper is to estimate what can be the effects of such resonances on the
orbital evolution of a satellite. For this, we evaluate precisely the capture
probabilities of a satellite in these resonances, and the variations in
orbital parameters due to their crossings in the case where the capture fails.
We also study the interaction with the evection resonance observed by Yokoyama
(2002) to determinate in which extent it influences the dynamics of a
satellite and determine the exact mechanism behind this excitation.

In Sect. 2, we first derive a secular model that is suitable to describe the
dynamics of a massless satellite of a rigid planet perturbed by the host star.
In Sect. 3, we identify all the possible secular resonances when the planet is
in a circular orbit around the star and its equator uniformly precesses with
constant obliquity. In Sect. 4, we study in detail with analytical models the
main eviction-like resonances that may significantly modify the inclination of
the satellite. In Sect. 5, we study the interaction of neighbor secular
resonances using frequency map analysis. In Sect. 6, we apply our model to the
future evolution of Phobos and numerically estimate the capture probabilities
in resonance for different values of the eccentricity and obliquity. Finally,
we discuss our results in Sect. 7.

##  2 Model

We consider a star orbited by a planet and a satellite, with masses $m_{0}$,
$m_{1}$ and $m$, respectively. The star and the satellite are considered point
masses. The planet is a rigid body with moments of inertia $A\leq B\leq C$,
and rotational angular momentum $\@vec{S}$.

The variables ($a$, $e$, $i$, $\lambda$, $\varpi$, $\Omega$) denote
respectively the semi-major axis, the eccentricity, the inclination with
respect to a reference plane, the mean longitude ($\lambda=\varpi+M$, with $M$
the mean anomaly), the longitude of the pericenter ($\varpi=\omega+\Omega$,
with $\omega$ the argument of the pericenter), and the longitude of the
ascending node of the orbit of the satellite around the planet. The same
variables with a subscript 1 are used for the orbit of the planet around the
star, that we name for simplicity, the ecliptic.

###  2.1 Secular Hamiltonian

Figure 1: Jacobi coordinates, where $\@vec{r}$ is the position of $m$ relative
to $m_{1}$ (satellite orbit), and $\@vec{r}_{1}$ is the position of the center
of mass of $m$ and $m_{1}$ relative to $m_{0}$ (planet orbit). The planet is a
rigid body, where $\@vec{S}$ is the rotational angular momentum.

We let $(\tilde{\@vec{R}}_{i},\@vec{R}_{i})$ be the barycentric canonical
variables with $\@vec{R}_{i}$ the position vector with respect to the
barycenter of the system, and $\tilde{\@vec{R}}_{i}=m_{i}\dot{\@vec{R}}_{i}$
the conjugate momenta ($i=0,1$, for the star, and planet, respectively, and we
do not use subscripts for the variables of the satellite). We perform a
canonical change into Jacobi variables $(\tilde{\@vec{r}}_{i},\@vec{r}_{i})$
as $\@vec{r}_{0}=\@vec{R}_{0}$,
$\@vec{r}_{1}=\delta\@vec{R}+(1-\delta)\@vec{R}_{1}-\@vec{R}_{0}$, and
$\@vec{r}=\@vec{R}-\@vec{R}_{1}$, with $\delta=m/(m_{1}+m)$ (see Fig. 1). In
the quadrupolar three-body problem approximation and with an average on the
fast angles of the rotation, the Hamiltonian of the system is given by (e.g.,
Smart 1953; Boué & Laskar 2006)

$\begin{split}\mathcal{H}=&\frac{\tilde{\@vec{r}}_{1}^{2}}{2\beta_{1}}+\frac{\tilde{\@vec{r}}^{2}}{2\beta}-\frac{\mathcal{G}m_{0}m_{1}}{r_{01}}-\frac{\mathcal{G}m_{0}m}{r^{\prime}}-\frac{\mathcal{G}m_{1}m}{r}\\\
&-\frac{\mathcal{E}m}{r^{3}}\left[1-3\left(\frac{\@vec{r}\cdot\@vec{s}}{r}\right)^{2}\right]-\frac{\mathcal{E}m_{0}}{r_{01}^{3}}\left[1-3\left(\frac{\@vec{r}_{01}\cdot\@vec{s}}{r_{01}}\right)^{2}\right]\
,\end{split}$ (1)

where $\mathcal{G}$ is the gravitational constant, $\beta=m_{1}m/(m_{1}+m)$,
$\beta_{1}=m_{0}(m_{1}+m)/(m_{0}+m_{1}+m)$,
$\@vec{r}_{01}=\@vec{R}_{1}-\@vec{R}_{0}=\@vec{r}_{1}-\delta\@vec{r}$,
$\@vec{r}^{\prime}=\@vec{R}-\@vec{R}_{0}=\@vec{r}_{1}+(1-\delta)\@vec{r}$,
$r=\|\@vec{r}\|$, $r^{\prime}=\|\@vec{r}^{\prime}\|$,
$r_{01}=\|\@vec{r}_{01}\|$, $\@vec{s}=\@vec{S}/\|\@vec{S}\|$, and

$\mathcal{E}=\frac{\mathcal{G}}{2}\left(C-\frac{A+B}{2}\right)\left(1-\frac{3}{2}\sin^{2}J\right)\
.$ (2)

The angle $J$ between the rotational angular momentum and the axis of maximal
inertia of the planet is very weak for the planets of the solar system, and so
we assume $\sin J=0$. We proceed to a development to the degree one in
$\delta$ and to the degree two in $r/r_{1}$. The Hamiltonian then becomes

$\begin{split}\mathcal{H}=&\frac{\tilde{\@vec{r}}_{1}^{2}}{2\beta_{1}}+\frac{\tilde{\@vec{r}}^{2}}{2\beta}-\frac{\mu_{1}\beta_{1}}{r_{1}}-\frac{\mu\beta}{r}+\frac{\mathcal{G}\beta
m_{0}}{2r_{1}^{3}}\left[r^{2}-\frac{3\left(\@vec{r}_{1}\cdot\@vec{r}\right)^{2}}{r_{1}^{2}}\right]\\\
&-\frac{\mathcal{E}m}{r^{3}}\left[1-3\left(\frac{\@vec{r}\cdot\@vec{s}}{r}\right)^{2}\right]-\frac{\mathcal{E}m_{0}}{r_{1}^{3}}\left[1-3\left(\frac{\@vec{r}_{1}\cdot\@vec{s}}{r_{1}}\right)^{2}\right]\\\
&-3\delta\frac{\mathcal{E}m_{0}}{r_{1}^{3}}\left[\frac{\@vec{r}_{1}\cdot\@vec{r}}{r_{1}^{2}}-5\frac{\left(\@vec{r}_{1}\cdot\@vec{r}\right)\left(\@vec{r}_{1}\cdot\@vec{s}\right)^{2}}{r_{1}^{4}}+2\frac{\left(\@vec{r}_{1}\cdot\@vec{s}\right)\left(\@vec{r}\cdot\@vec{s}\right)}{r_{1}^{2}}\right]\
,\end{split}$ (3)

with $\mu=\mathcal{G}(m_{1}+m)$ and $\mu_{1}=\mathcal{G}(m_{0}+m_{1}+m)$.

We now consider that the satellite is a massless particle. We thus perform a
canonical change of variables
$(\tilde{\@vec{r}},\@vec{r})\rightarrow(\tilde{\@vec{r}}^{\prime},\@vec{r})$,
with $\tilde{\@vec{r}}^{\prime}=\tilde{\@vec{r}}/m$, and then make
$m\rightarrow 0$, in order to obtain the Hamiltonian for the motion of the
perturbed satellite

$\begin{split}\mathcal{H}=&\frac{\tilde{\@vec{r}}^{\prime
2}}{2}-\frac{\mu}{r}+\frac{\mathcal{G}m_{0}}{2r_{1}^{3}}\left[r^{2}-\frac{3\left(\@vec{r}_{1}\cdot\@vec{r}\right)^{2}}{r_{1}^{2}}\right]-\frac{\mathcal{E}}{r^{3}}\left[1-3\left(\frac{\@vec{r}\cdot\@vec{s}}{r}\right)^{2}\right]\\\
&-\frac{3\mathcal{E}}{r_{1}^{3}}\frac{m_{0}}{m_{1}}\left[\frac{\@vec{r}_{1}\cdot\@vec{r}}{r_{1}^{2}}-5\frac{\left(\@vec{r}_{1}\cdot\@vec{r}\right)\left(\@vec{r}_{1}\cdot\@vec{s}\right)^{2}}{r_{1}^{4}}+2\frac{\left(\@vec{r}_{1}\cdot\@vec{s}\right)\left(\@vec{r}\cdot\@vec{s}\right)}{r_{1}^{2}}\right].\end{split}$
(4)

Finally, we average the orbit of the satellite over the mean anomaly (for
details see Boué & Laskar 2006), and obtain the secular Hamiltonian

$\begin{split}\mathcal{H}_{s}=&\frac{3\mathcal{G}m_{0}a^{2}}{2r_{1}^{3}}\left[\@vec{e}^{2}+\frac{1}{2}\frac{\left(\@vec{r}_{1}\cdot\@vec{u}\right)^{2}}{r_{1}^{2}}-\frac{5}{2}\frac{\left(\@vec{r}_{1}\cdot\@vec{e}\right)^{2}}{r_{1}^{2}}\right]\\\
&+\frac{\mathcal{E}}{2a^{3}\|\@vec{u}\|^{5}}\left[\@vec{u}^{2}-3\left(\@vec{s}\cdot\@vec{u}\right)^{2}\right]\\\
&+\frac{9a\mathcal{E}}{r_{1}^{3}}\frac{m_{0}}{m_{1}}\left[\frac{\@vec{r}_{1}\cdot\@vec{e}}{2r_{1}^{2}}\left(1-5\frac{\left(\@vec{r}_{1}\cdot\@vec{s}\right)^{2}}{r_{1}^{2}}\right)+\frac{\left(\@vec{r}_{1}\cdot\@vec{s}\right)\left(\@vec{e}\cdot\@vec{s}\right)}{r_{1}^{2}}\right]\
,\end{split}$ (5)

with

$\@vec{e}=e\,\@vec{i}\
,\quad\mathrm{and}\quad\@vec{u}=\sqrt{1-e^{2}}\,\@vec{k}\ ,$ (6)

where $\@vec{k}=\@vec{G}/G$ is the normal to the orbit of the satellite,
$\@vec{i}$ is the unit vector indicating the direction of its pericenter,
$\@vec{G}$ is the orbital angular momentum, $G=L\sqrt{1-e^{2}}$, and
$L=\sqrt{\mathcal{G}m_{1}a}$. The first term of $\mathcal{H}_{s}$ (Eq. (5))
corresponds to the perturbations from the central star, the second term
corresponds to perturbations from the gravitational flattening of the planet,
while the third term describes the effect of the couple exerted by the central
star on the rigid planet on the orbital motion of the satellite.

Figure 2: Orbital elements of the orbit of the satellite (in red).

###  2.2 Equations of motion

The equations of motion for the vectors $\@vec{e}$ and $\@vec{u}$ can be
obtained from the secular Hamiltonian $\mathcal{H}_{s}$ (Eq. (5)) as (Tremaine
et al. 2009; Farago & Laskar 2010)

$\begin{cases}\dfrac{d\@vec{e}}{dt}=\dfrac{1}{L}\left(\nabla_{\@vec{u}}\mathcal{H}_{s}\times\@vec{e}+\nabla_{\@vec{e}}\mathcal{H}_{s}\times\@vec{u}\right)\
,\cr\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\dfrac{d\@vec{u}}{dt}=\dfrac{1}{L}\left(\nabla_{\@vec{u}}\mathcal{H}_{s}\times\@vec{u}+\nabla_{\@vec{e}}\mathcal{H}_{s}\times\@vec{e}\right)\
,\end{cases}$ (7)

which gives

$\begin{cases}\dfrac{d\@vec{e}}{dt}=\@vec{h}_{1}\times\@vec{e}+\@vec{h}_{2}\times\@vec{u}+\gamma\,\@vec{u}\times\@vec{e}\
,\cr\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\dfrac{d\@vec{u}}{dt}=\@vec{h}_{1}\times\@vec{u}+\@vec{h}_{2}\times\@vec{e}\
,\end{cases}$ (8)

with

$\begin{split}\@vec{h}_{1}=&\frac{1}{L}\left[\frac{3\mathcal{G}m_{0}a^{2}}{2r_{1}^{5}}\left(\@vec{r}_{1}\cdot\@vec{u}\right)\@vec{r}_{1}-\frac{3\mathcal{E}}{a^{3}\|\@vec{u}\|^{5}}\left(\@vec{s}\cdot\@vec{u}\right)\@vec{\@vec{s}}\right]\
,\end{split}$
$\begin{split}\@vec{h}_{2}=&\frac{1}{L}\left[-\frac{15\mathcal{G}m_{0}a^{2}}{2r_{1}^{5}}\left(\@vec{r}_{1}\cdot\<EMAIL_ADDRESS>&\left.+\frac{9\mathcal{E}am_{0}}{r_{1}^{5}m_{1}}\left(\frac{\@vec{r}_{1}}{2}-\frac{5}{2}\frac{\left(\@vec{r}_{1}\cdot\@vec{s}\right)^{2}}{r_{1}^{2}}\@vec{\@vec{r}_{1}}+\left(\@vec{r}_{1}\cdot\@vec{s}\right)\@vec{s}\right)\right]\
,\end{split}$
$\gamma=-\frac{1}{L}\left[\frac{3\mathcal{E}}{2a^{3}\|\@vec{u}\|^{5}}\left(1-5\left(\@vec{s}\cdot\@vec{k}\right)^{2}\right)+\frac{3\mathcal{G}m_{0}a^{2}}{r_{1}^{3}}\right]\
.$

We use Eq. (8) to obtain numerically the secular orbital evolution of the
satellite. To completely solve these equations, we also need to compute the
motion of the planet around the star, $\@vec{r}_{1}(t)$, and the evolution of
the spin axis, $\@vec{s}(t)$. For simplicity, we assume that the planet moves
on a circular orbit and that most of the angular momentum is on its orbit. We
also assume that the spin of the planet can only precess around the normal to
the ecliptic with constant speed and obliquity (i.e., the angle between the
equator of the planet and the ecliptic). Therefore, we have

$\@vec{r}_{1}\left(t\right)=\mathcal{R}_{3}\left(n_{1}t\right)\begin{pmatrix}a_{1}\\\
0\\\ 0\\\ \end{pmatrix}\ ,$ (9)

and

$\@vec{s}\left(t\right)=\mathcal{R}_{3}\left(\dot{\psi}t+\psi\left(0\right)\right)\mathcal{R}_{1}\left(\varepsilon\right)\begin{pmatrix}0\\\
0\\\ 1\\\ \end{pmatrix}\ ,$ (10)

where $n_{1}$ is the mean motion, $\varepsilon$ is the obliquity, and $\psi$
is the precession angle, which measures the precession of the equator with
respect to the ecliptic (Fig. 2). $\mathcal{R}_{1}$ and $\mathcal{R}_{3}$ are
the rotations about the $x-$axis and $z-$axis, respectively

$\mathcal{R}_{1}\left(\alpha\right)=\begin{pmatrix}1&0&0\\\
0&\cos\alpha&-\sin\alpha\\\ 0&\sin\alpha&\cos\alpha\\\ \end{pmatrix}\ ,$ (11)
$\mathcal{R}_{3}\left(\alpha\right)=\begin{pmatrix}\cos\alpha&-\sin\alpha&0\\\
\sin\alpha&\cos\alpha&0\\\ 0&0&1\\\ \end{pmatrix}\ .$ (12)

##  3 Secular resonances

In general, the last term of the secular Hamiltonian (Eq. (5)) is much weaker
than the other two terms. For instance, in the case of the Martian satellite
Phobos, the ratio between the amplitudes of the third and first term is $\sim
10^{-9}$. Moreover, on average, the effect of this third term is null over a
satellite pericenter precession period. Therefore, we can often neglect its
contribution, and the secular Hamiltonian (Eq. (5)) is rewritten as

$\begin{split}\mathcal{H}=&\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left[e^{2}+\frac{1-e^{2}}{2}\left(\@vec{\hat{r}}_{1}\cdot\@vec{k}\right)^{2}-\frac{5e^{2}}{2}\left(\@vec{\hat{r}}_{1}\cdot\@vec{i}\right)^{2}\right]\\\
&+\frac{\mathcal{E}}{2a^{3}\left(1-e^{2}\right)^{3/2}}\left[1-3\left(\@vec{s}\cdot\@vec{k}\right)^{2}\right]\
,\end{split}$ (13)

with $\@vec{\hat{r}}_{1}=\@vec{r}_{1}/r_{1}$. Resonances between the orbit of
the satellite and the Sun can only occur due to the first term. In order to
identify all possible resonances, we need to develop $\mathcal{H}$ in
trigonometric series. This development is given by Yokoyama (2002) but with
some misprints. A corrected version is presented in Carvalho & Vilhena de
Moraes (2020).

We consider here the ecliptic as reference plane (Fig. 2). Assuming that the
satellite is a massless particle and that the orbital angular momentum of the
planet is much larger than its rotational angular momentum, the ecliptic
remains an inertial plane to a very good approximation. We then need to
express $(\@vec{\hat{r}}_{1}\cdot\@vec{k})^{2}$ and
$(\@vec{\hat{r}}_{1}\cdot\@vec{i})^{2}$ as a function of the orbital elements
of the satellite. The unit vectors $\@vec{k}$, $\@vec{i}$ and
$\@vec{\hat{r}}_{1}$ can be expressed with respect to the ecliptic as a
succession of rotations

$\@vec{k}=\mathcal{R}_{3}\left(\psi\right)\mathcal{R}_{1}\left(\varepsilon\right)\mathcal{R}_{3}\left(\Omega\right)\mathcal{R}_{1}\left(i\right)\begin{pmatrix}0\\\
0\\\ 1\\\ \end{pmatrix}\ ,$ (14)
$\@vec{i}=\mathcal{R}_{3}\left(\psi\right)\mathcal{R}_{1}\left(\varepsilon\right)\mathcal{R}_{3}\left(\Omega\right)\mathcal{R}_{1}\left(i\right)\mathcal{R}_{3}\left(\omega\right)\begin{pmatrix}1\\\
0\\\ 0\\\ \end{pmatrix}\ ,$ (15)

and

$\@vec{\hat{r}}_{1}=\mathcal{R}_{3}\left(\lambda_{1}\right)\begin{pmatrix}1\\\
0\\\ 0\\\ \end{pmatrix}\ ,$ (16)

where $\lambda_{1}=n_{1}t$ is the mean longitude of the planet (we assume a
circular orbit). The Hamiltonian (Eq. (13)) then becomes

$\begin{split}&\mathcal{H}=\frac{\mathcal{E}}{2a^{3}\left(1-e^{2}\right)^{3/2}}\left[1-3\cos^{2}i\right]+\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left[-\frac{e^{2}}{4}\right.\\\
&+\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\left(\sin^{2}i+\sin^{2}\varepsilon-\frac{3}{2}\sin^{2}i\sin^{2}\varepsilon\right)\\\
&-\frac{1}{8}\left(1+\frac{3}{2}e^{2}\right)\sin^{2}i\sin^{2}\varepsilon\cos\left(2\Omega\right)\\\
&+\frac{1}{8}\left(1+\frac{3}{2}e^{2}\right)\sin\left(2i\right)\sin\left(2\varepsilon\right)\cos\left(\Omega\right)\\\
&-\frac{5}{8}e^{2}\sin^{2}i\left(1-\frac{3}{2}\sin^{2}\varepsilon\right)\cos\left(2\omega\right)\\\
&-\frac{5}{32}e^{2}\left(1+\cos
i\right)^{2}\sin^{2}\varepsilon\cos\left(2\left(\Omega+\omega\right)\right)\\\
&-\frac{5}{32}e^{2}\left(1-\cos
i\right)^{2}\sin^{2}\varepsilon\cos\left(2\left(\Omega-\omega\right)\right)\\\
&-\frac{5}{16}e^{2}\sin i\left(1+\cos
i\right)\sin\left(2\varepsilon\right)\cos\left(\Omega+2\omega\right)\\\
&+\frac{5}{16}e^{2}\sin i\left(1-\cos
i\right)\sin\left(2\varepsilon\right)\cos\left(\Omega-2\omega\right)\\\
&-\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\left(1-\frac{3}{2}\sin^{2}i\right)\sin^{2}\varepsilon\cos\left(2\left(\lambda_{1}-\psi\right)\right)\\\
&-\frac{1}{16}\left(1+\frac{3}{2}e^{2}\right)\sin^{2}i\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+\Omega\right)\right)\\\
&-\frac{1}{16}\left(1+\frac{3}{2}e^{2}\right)\sin^{2}i\left(1+\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi-\Omega\right)\right)\\\
&+\frac{1}{8}\left(1+\frac{3}{2}e^{2}\right)\sin\left(2i\right)\sin\varepsilon\left(1-\cos\varepsilon\right)\cos\left(2\left(\lambda_{1}-\psi\right)+\Omega\right)\\\
&-\frac{1}{8}\left(1+\frac{3}{2}e^{2}\right)\sin\left(2i\right)\sin\varepsilon\left(1+\cos\varepsilon\right)\cos\left(2\left(\lambda_{1}-\psi\right)-\Omega\right)\\\
&-\frac{15}{32}e^{2}\sin^{2}i\sin^{2}\varepsilon\cos\left(2\left(\lambda_{1}-\psi+\omega\right)\right)\\\
&-\frac{15}{32}e^{2}\sin^{2}i\sin^{2}\varepsilon\cos\left(2\left(\lambda_{1}-\psi-\omega\right)\right)\\\
&-\frac{5}{64}e^{2}\left(1+\cos
i\right)^{2}\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+\Omega+\omega\right)\right)\\\
&-\frac{5}{64}e^{2}\left(1-\cos
i\right)^{2}\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+\Omega-\omega\right)\right)\\\
&-\frac{5}{64}e^{2}\left(1-\cos
i\right)^{2}\left(1+\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi-\Omega+\omega\right)\right)\\\
&-\frac{5}{64}e^{2}\left(1+\cos
i\right)^{2}\left(1+\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi-\Omega-\omega\right)\right)\\\
&-\frac{5}{16}e^{2}\sin i\left(1+\cos
i\right)\sin\varepsilon\left(1-\cos\varepsilon\right)\\\
&\cos\left(2\left(\lambda_{1}-\psi\right)+\Omega+2\omega\right)\\\
&+\frac{5}{16}e^{2}\sin i\left(1-\cos
i\right)\sin\varepsilon\left(1-\cos\varepsilon\right)\\\
&\cos\left(2\left(\lambda_{1}-\psi\right)+\Omega-2\omega\right)\\\
&-\frac{5}{16}e^{2}\sin i\left(1-\cos
i\right)\sin\varepsilon\left(1+\cos\varepsilon\right)\\\
&\cos\left(2\left(\lambda_{1}-\psi\right)-\Omega+2\omega\right)\\\
&+\frac{5}{16}e^{2}\sin i\left(1+\cos
i\right)\sin\varepsilon\left(1+\cos\varepsilon\right)\\\
&\left.\cos\left(2\left(\lambda_{1}-\psi\right)-\Omega-2\omega\right)\right].\end{split}$
(17)

From this Hamiltonian, we can identify 14 different possible combinations of
angles between the orbital motion of the planet and the secular motion of the
satellite. The corresponding possible resonances are then

$\displaystyle\nu_{e}\equiv n^{\prime}_{1}-\dot{\Omega}-\dot{\omega}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
n^{\prime}_{1}+\dot{\Omega}+\dot{\omega}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle n^{\prime}_{1}-\dot{\Omega}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\nu_{1}\equiv n^{\prime}_{1}+\dot{\Omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle n^{\prime}_{1}-\dot{\omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle n^{\prime}_{1}+\dot{\omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle n^{\prime}_{1}+\dot{\Omega}-\dot{\omega}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
n^{\prime}_{1}-\dot{\Omega}+\dot{\omega}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle 2n^{\prime}_{1}-\dot{\Omega}$ $\displaystyle=$ $\displaystyle
0$ $\displaystyle\nu_{2}\equiv 2n^{\prime}_{1}+\dot{\Omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle 2n^{\prime}_{1}+\dot{\Omega}+2\dot{\omega}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
2n^{\prime}_{1}+\dot{\Omega}-2\dot{\omega}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle 2n^{\prime}_{1}-\dot{\Omega}+2\dot{\omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\nu_{i}\equiv
2n^{\prime}_{1}-\dot{\Omega}-2\dot{\omega}$ $\displaystyle=$ $\displaystyle 0\
,$ (18)

where $n^{\prime}_{1}=n_{1}-\dot{\psi}$ is the mean motion of the planet
corrected by the precession frequency of the spin axis.

The argument of the pericenter, $\omega$, is not an inertial angle, as it is
defined with respect to the line of nodes between the orbit of the satellite
and the equator of the planet (Fig. 2), but it can be related with the
inertial longitudes $\omega=\varpi-\Omega$. The first resonance in the list
(18) thus becomes $\nu_{e}=n^{\prime}_{1}-\dot{\varpi}=0$, and occurs when the
pericenter precession frequency of the satellite, $\dot{\varpi}$, is equal to
$n^{\prime}_{1}$. This corresponds to the well-known evection resonance. When
capture in this resonance occurs, we can observe a significant increase in the
eccentricity of the satellite (e.g., Touma & Wisdom 1998; Frouard et al.
2010). However, the amplitude of the evection resonance is proportional to
$e^{2}$, and so this resonance cannot occur for satellites in nearly circular
orbits, which is the final outcome of tidal evolution. The resonance
$n^{\prime}_{1}-\dot{\omega}=0$ has been studied in the case of Phobos by
Yoder (1982), who showed that its crossing could lead to a sudden increase in
the eccentricity of Phobos of maximal amplitude $0.016$. We also note that in
expression (17) other secular resonances not involving the mean motion $n_{1}$
are also possible, such as $\dot{\omega}=0$ (Kozai resonance), but we do not
study them here.

##  4 Eviction-like resonances

The last resonance in the list (18),
$\nu_{i}=2n^{\prime}_{1}-2\dot{\varpi}+\dot{\Omega}=0$, corresponds to the
”eviction” resonance described by Touma & Wisdom (1998), who showed that it
can increase the inclination of the satellite. In principle, any resonance
containing the longitude of the ascending node, $\Omega$, can excite the
inclination, and hence be dubbed as ”eviction-like” resonance. However, as for
the evection resonance, the amplitudes of most of these resonances are also
proportional to the square of the eccentricity (Eq. (17)), and so they are not
present for nearly circular satellite orbits. This incidentally includes the
”original” eviction resonance $\nu_{i}$ reported by Touma & Wisdom (1998).

In the Hamiltonian (17), only four resonances,
$n^{\prime}_{1}\pm\dot{\Omega}=0$ and $2n^{\prime}_{1}\pm\dot{\Omega}=0$, do
not depend on the pericenter precession frequency, and thus do not have a null
amplitude for an eccentricity equal to zero. As a consequence, they are able
to modify the inclination of a satellite even if the eccentricity is nearly
zero, as currently observed for all tidally evolved satellites in the solar
system. Nonetheless, since $n_{1}>0$ and $\dot{\Omega}$ is expected to be
negative for a prograde orbit, only the resonances
$\nu_{1}=n^{\prime}_{1}+\dot{\Omega}=0$ and
$\nu_{2}=2n^{\prime}_{1}+\dot{\Omega}=0$ can occur (Yokoyama 2002). Therefore,
in this section, we focus our study on these two resonances.

In general, we have $|\dot{\psi}|\ll n_{1}$, thus, for simplicity, we assume
$\dot{\psi}=0$ only in this section, such that we can use the Delaunay
canonical variables. It is possible to obtain these variables for a moving
equator, but we needed to modify the Hamiltonian (e.g., Goldreich 1965;
Kinoshita 1993).

###  4.1 Resonance $\nu_{1}$

We denote $\nu_{1}$ the resonance associated to the relation
$n^{\prime}_{1}+\dot{\Omega}=0$. Among the resonant terms of the Hamiltonian
(17), we keep only the term associated with $\nu_{1}$, and we get for the
resonant Hamiltonian

$\begin{split}\mathcal{H}_{\nu_{1}}=&\frac{\mathcal{E}}{2a^{3}\left(1-e^{2}\right)^{3/2}}\left(1-3\cos^{2}i\right)+\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left(-\frac{e^{2}}{4}\right.\\\
&+\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\left(\sin^{2}i+\sin^{2}\varepsilon-\frac{3}{2}\sin^{2}i\sin^{2}\varepsilon\right)\\\
&\left.-\frac{1}{16}\left(1+\frac{3}{2}e^{2}\right)\sin^{2}i\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+\Omega\right)\right)\right).\end{split}$
(19)

Using the Delaunay variables for the satellite $(L,G,H,l,g,h)$, where
$L=\sqrt{\mathcal{G}m_{1}a}$, $G=L\sqrt{1-e^{2}}$, $H=G\cos i$, $l=M$,
$g=\omega$, and $h=\Omega$, we obtain

$\begin{split}&\mathcal{H}_{\nu_{1}}=\frac{\mathcal{E}\left(\mathcal{G}m_{1}\right)^{3/2}}{2a^{3/2}G^{3}}\left(1-3\frac{H^{2}}{G^{2}}\right)+\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left(-\frac{L^{2}-G^{2}}{4L^{2}}\right.\\\
&+\frac{1}{4}\left(1+\frac{3}{2}\frac{L^{2}-G^{2}}{L^{2}}\right)\left(\sin^{2}\varepsilon+\left(1-\frac{H^{2}}{G^{2}}\right)\left(1-\frac{3}{2}\sin^{2}\varepsilon\right)\right)\\\
&-\frac{1}{16}\left(1+\frac{3}{2}\frac{L^{2}-G^{2}}{L^{2}}\right)\left(1-\frac{H^{2}}{G^{2}}\right)\\\
&\left.\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+h\right)\right)\right).\end{split}$
(20)

The Hamiltonian $\mathcal{H}_{\nu_{1}}$ does not depend on the angles $l$ and
$g$, and so the associated action variables $L$ and $G$ are constant. As a
result, the semi-major axis and the eccentricity are also constant, that is,
the $\nu_{1}$ resonance does not modify the eccentricity. By removing the
constant terms and the terms which do not depend on $H$ or $h$, the
Hamiltonian becomes

$\mathcal{H}_{\nu_{1}}=-\mathcal{A}H^{2}-\mathcal{B}\left(G^{2}-H^{2}\right)\cos\left(2\left(\lambda_{1}-\psi+h\right)\right),$
(21)

with

$\mathcal{A}=\frac{3\mathcal{G}m_{0}a^{2}}{8G^{2}a_{1}^{3}}\left(1+\frac{3}{2}e^{2}\right)\left(1-\frac{3}{2}\sin^{2}\varepsilon\right)+\frac{3\mathcal{E}}{2G^{2}a^{3}\left(1-e^{2}\right)^{3/2}},$
(22)
$\mathcal{B}=\frac{3\mathcal{G}m_{0}a^{2}}{32G^{2}a_{1}^{3}}\left(1+\frac{3}{2}e^{2}\right)\left(1-\cos\varepsilon\right)^{2}.$
(23)

As in Touma & Wisdom (1998) and Yokoyama (2002), we use a generating function
to perform a canonical change of variables, $F_{1}=(\lambda_{1}-\psi+h)X$. We
obtain the canonical variables $(X,x)$ with $X=H$ and $x=\lambda_{1}-\psi+h$.
As we suppose $\dot{\psi}=0$ in this section, we have
$\partial(\lambda_{1}-\psi)/\partial t=n_{1}$, and the Hamiltonian becomes

$\mathcal{H}_{\nu_{1}}=-\mathcal{A}X^{2}-\mathcal{B}\left(G^{2}-X^{2}\right)\cos
2x+n_{1}X\ .$ (24)

The equations of motion are then given by

$\dot{X}=-\frac{\partial\mathcal{H}_{\nu_{1}}}{\partial
x}=-2\mathcal{B}\left(G^{2}-X^{2}\right)\sin 2x\ ,$ (25)
$\dot{x}=\frac{\partial\mathcal{H}_{\nu_{1}}}{\partial
X}=-2\mathcal{A}X+2\mathcal{B}X\cos 2x+n_{1}\ .$ (26)

This Hamiltonian possesses four fixed points, obtained when
$(\dot{X},\dot{x})=(0,0)$: the point
$(X=n_{1}/(2(\mathcal{A}+\mathcal{B})),x=\pi/2)$ and the point
$(X=n_{1}/(2(\mathcal{A}+\mathcal{B})),x=3\pi/2)$ are stable, while the points
$(X=n_{1}/(2(\mathcal{A}-\mathcal{B})),x=0)$ and
$(X=n_{1}/(2(\mathcal{A}-\mathcal{B})),x=\pi)$ are unstable. In Fig. 3 we show
the level curves of the Hamiltonian (24) for Phobos, assuming an obliquity
$\varepsilon=90^{\circ}$ for Mars (Table 1), where we can clearly identify the
fixed points as well as the separatrix that encircles the resonant area. If
the satellite is at a stable fixed point, its inclination verifies

$\cos i=\frac{n_{1}}{2G\left(\mathcal{A}+\mathcal{B}\right)}\ ,$ (27)

while for an unstable fixed point, the inclination is given by

$\cos i=\frac{n_{1}}{2G\left(\mathcal{A}-\mathcal{B}\right)}\ .$ (28)

Therefore, the stable points are present only if $n_{1}\leq
2G(\mathcal{A}+\mathcal{B})$, and the unstable ones if $n_{1}\leq
2G(\mathcal{A}-\mathcal{B})$. The two equilibrium points are usually very
close to each other. Indeed, near the $\nu_{1}$ resonance, the perturbation
due to the flattening of the planet is in general much larger than the stellar
perturbation, and so

$\mathcal{A}+\mathcal{B}\approx\mathcal{A}-\mathcal{B}\approx\frac{3\mathcal{E}}{2G^{2}a^{3}(1-e^{2})^{3/2}}\
.$ (29)

In Fig. 4, we show the stable equilibria for Phobos computed with the values
from Table 1, assuming an obliquity of Mars similar to the present one,
$\varepsilon=25^{\circ}$. We observe that there is a critical value for the
semi-major axis of the satellite (obtained with $\cos i=1$),

$a_{\nu_{1}}\approx\left(\frac{3\mathcal{E}}{\sqrt{\mathcal{G}m_{1}}n_{1}(1-e^{2})^{2}}\right)^{2/7},$
(30)

above which the $\nu_{1}$ resonance is not present. In the case of Phobos, we
have $a_{\nu_{1}}/R\approx 2.617$. This value is very close to the surface of
Mars, but still larger than the Roche limit for a solid body with the density
of Phobos, $a/R\approx 1.61$ (e.g., Chandrasekhar 1987).

Figure 3: Level curves of the Hamiltonian (24) for Phobos (Table 1), for a
semi-major axis $a=2.61686\,R$, and assuming an obliquity
$\varepsilon=90^{\circ}$ for Mars. Table 1: Values of the present parameters
of Phobos, Mars, and the Sun used in this work. For some parameters, we give
an approximate value.

parameter | unit | value
---|---|---
$\mathcal{G}m_{0}$ | $\mathrm{au}^{3}/\mathrm{day}^{2}$ | $2.959122\times 10^{-4}$
$\mathcal{G}m_{1}$ | $\mathrm{au}^{3}/\mathrm{day}^{2}$ | $9.5495351\times 10^{-11}$
$\mathcal{G}m\,^{(a)}$ | $\mathrm{km}^{3}/\mathrm{s}^{2}$ | $7.092\times 10^{-4}$
$a$ | $\mathrm{km}$ | $9375.$
$e$ | $-$ | $0.015$
$i$ | $\mathrm{\SIUnitSymbolDegree}$ | $1.$
$T_{1}$ | $\mathrm{day}$ | $686.98$
$n_{1}$ | rad/day | $2\pi/T_{1}$
$a_{1}$ | $\mathrm{au}$ | $(\mathcal{G}(m_{0}+m_{1})/n_{1}^{2})^{1/3}$
$\varepsilon$ | $\mathrm{\SIUnitSymbolDegree}$ | $25.$
$\dot{\psi}\,^{(b)}$ | $\mathrm{\SIUnitSymbolArcsecond}\mathrm{/}\mathrm{yr}$ | $-7.6083$
$J_{2}\,^{(b)}$ | $-$ | $1.95661\times 10^{-3}$
$R\,^{(b)}$ | $\mathrm{km}$ | $3396.$
$k_{2}\,^{(c)}$ | $-$ | $0.183$
$Q\,^{(c)}$ | $-$ | $99.5$
111$J_{2}=(C-\frac{1}{2}(A+B))/(m_{1}R^{2})$; (a) Jacobson (2010); (b)
Konopliv et al. (2016); (c) Jacobson & Lainey (2014).

Figure 4: Evolution of the inclination of Phobos for the stable fixed points
with respect to its semi-major axis for the $\nu_{1}$ resonance (red curve)
computed with Eq. (27), and for the $\nu_{2}$ resonance (blue curve) computed
with Eq. (37) for a Martian obliquity of
$\varepsilon=25\mathrm{\SIUnitSymbolDegree}$. The Roche limit for Phobos is
$a/R\approx 1.6$.

###  4.2 Resonance $\nu_{2}$

We denote $\nu_{2}$ the resonance associated to the relation
$2n^{\prime}_{1}+\dot{\Omega}=0$. Among the resonant terms of the Hamiltonian
(17), we keep only the term associated with $\nu_{2}$, and we get for the
resonant Hamiltonian

$\begin{split}&\mathcal{H}_{\nu_{2}}=\frac{\mathcal{E}}{2a^{3}\left(1-e^{2}\right)^{3/2}}\left(1-3\cos^{2}i\right)+\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left(-\frac{e^{2}}{4}\right.\\\
&+\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\left(\sin^{2}i+\sin^{2}\varepsilon-\frac{3}{2}\sin^{2}i\sin^{2}\varepsilon\right)\\\
&\left.+\frac{1}{8}\left(1+\frac{3}{2}e^{2}\right)\sin\left(2i\right)\sin\varepsilon\left(1-\cos\varepsilon\right)\cos\left(2\left(\lambda_{1}-\psi\right)+\Omega\right)\right).\end{split}$
(31)

As for the $\nu_{1}$ resonance, the Hamiltonian does not depend on the angle
$g$ and the eccentricity remains constant. The Hamiltonian then becomes

$\mathcal{H}_{\nu_{2}}=-\mathcal{A}H^{2}+\mathcal{C}H\sqrt{G^{2}-H^{2}}\cos\left(2\left(\lambda_{1}-\psi\right)+h\right),$
(32)

where $\mathcal{A}$ is given by expression (22), and

$\mathcal{C}=\frac{3\mathcal{G}m_{0}a^{2}}{8G^{2}a_{1}^{3}}\left(1+\frac{3}{2}e^{2}\right)\sin\varepsilon\left(1-\cos\varepsilon\right).$
(33)

We now consider the generating function $F_{2}=(2(\lambda_{1}-\psi)+h)Y$. The
canonical variables $(Y,y)$ verify $Y=H$ and $y=2(\lambda_{1}-\psi)+h$, and
the resonant Hamiltonian becomes

$\mathcal{H}_{\nu_{2}}=-\mathcal{A}Y^{2}+\mathcal{C}Y\sqrt{G^{2}-Y^{2}}\cos
y+2n_{1}Y\ .$ (34)

The equations of motion are then given by

$\dot{Y}=-\frac{\partial\mathcal{H}_{\nu_{2}}}{\partial
y}=\mathcal{C}Y\sqrt{G^{2}-Y^{2}}\sin y,$ (35)
$\dot{y}=\frac{\partial\mathcal{H}_{\nu_{2}}}{\partial
Y}=-2\mathcal{A}Y+\mathcal{C}\frac{G^{2}-2Y^{2}}{\sqrt{G^{2}-Y^{2}}}\cos
y+2n_{1}\ .$ (36)

The fixed points are then located at $y=0$ or $y=\pi$. In order to obtain the
$Y$ values for the fixed points, we must solve the following quartic equation

$z^{4}-2bcz^{3}-\left(1-b^{2}\right)z^{2}+2bcz+\frac{1-c^{2}}{4}-b^{2}=0,$
(37)

with $z=Y/G$, $c=\mathcal{A}/\sqrt{\mathcal{A}^{2}+\mathcal{C}^{2}}$, and
$b=n_{1}/(G\sqrt{\mathcal{A}^{2}+\mathcal{C}^{2}})$. Contrary to the $\nu_{1}$
resonance, it is not possible to find a simple analytical expression for the
fixed points, although we can find the roots numerically. In Fig. 4, we show
the stable fixed point located at $y=\pi$ for Phobos with the values from
Table 1, and assuming an obliquity of Mars similar to the present one,
$\varepsilon=25^{\circ}$. For a satellite close to the $\nu_{2}$ resonance,
the perturbations due to the flattening of the planet are in general larger
than the stellar perturbations, and so $\mathcal{A}\approx
3\mathcal{E}/(2G^{2}a^{3}(1-e^{2})^{3/2})\gg\mathcal{C}$. It is then possible
to obtain an approximate value for the inclination of the fixed points. By
neglecting the term proportional to $\mathcal{C}$ in Eq. (36), we obtain

$-2\mathcal{A}Y+2n_{1}=0\ ,$ (38)

and the inclination is given by

$\cos i\approx\frac{n_{1}}{G\mathcal{A}}\ .$ (39)

Therefore, as for the $\nu_{1}$ resonance, there is a critical value
$a_{\nu_{2}}$ of the semi-major axis above which the $\nu_{2}$ resonance is
not present. An approximate value of $a_{\nu_{2}}$ is given by (obtained with
$\cos i=1$)

$a_{\nu_{2}}\approx\left(\frac{3\mathcal{E}}{2\sqrt{\mathcal{G}m_{1}}n_{1}(1-e^{2})^{2}}\right)^{2/7}\approx\frac{a_{\nu_{1}}}{2^{2/7}}\approx
0.82\,a_{\nu_{1}}\ ,$ (40)

where $a_{\nu_{1}}$ is the critical value for the $\nu_{1}$ resonance (Eq.
(30)). In the case of Phobos, we have $a_{\nu_{2}}/R\approx 2.147$.

###  4.3 Capture probabilities

Due to tidal effects, the semi-major axis of the satellites usually evolves
with time. When the semi-major axis is increasing, the eviction-like
resonances are crossed if we initially have $a<a_{\nu}$. As a result, capture
is not possible, because the resonances are no longer present for $a>a_{\nu}$
(Fig. 4). As noted by Touma & Wisdom (1998), the crossing of the resonance can
nevertheless lead to a sudden variation in the inclination of the satellite
due to the conservation of the area in the phase plane (Fig. 3).

On the other hand, when the semi-major axis is decreasing, the eviction-like
resonances are encountered if we initially have $a>a_{\nu}$. As a result,
capture in resonance may occur. In that case, as the semi-major axis is
decreasing, the inclination increases (Fig. 4). We then conclude that, for a
satellite with an orbit close to the equator, capture in eviction-like
resonances is only possible if the semi-major axis is decreasing.

It is possible to compute analytically the capture probability in the
$\nu_{1}$ resonance, provided that the evolution of the semi-major axis is
adiabatic. To be considered as adiabatic, the evolution of the semi-major axis
must be much slower than the conservative inclination variations, that is,
tidal dissipation must be weak. In general, in the phase space of a given
resonance we can distinguish three different regions delimited by the
separatrix of the resonance (Fig. 3): a libration zone, with area $A_{lib}$,
and two circulation zones, one above the libration zone, with area $A_{circ}$,
and another below the libration zone, with area $A_{circ}^{\prime}$. The
capture probability is then obtained by the modification of the phase space
with time, that is, by the change in the areas encircled by the separatrix
(Yoder 1979; Henrard 1982, 1993)

$P_{\mathrm{cap}}=\dot{A}_{lib}/\dot{A}_{tot}\ ,$ (41)

where $A_{tot}$ corresponds to the sum of the libration area with the
circulation area increasing with time. The ordinate of the stable fixed point
inside the libration area is (Eq. (27))

$X=\frac{n_{1}}{2(\mathcal{A}+\mathcal{B})}\approx\frac{n_{1}G^{2}a^{3}(1-e^{2})^{3/2}}{3\mathcal{E}}\propto
a^{4}\ ,$ (42)

and so the ordinate of the fixed point decreases when the semi-major
decreases, which means that the circulation area whose surface increases with
time is $A_{circ}$, and then $A_{tot}=A_{lib}+A_{circ}$.

As the phase space of the $\nu_{1}$ resonance is $\pi-$periodic (Eq. (24)), we
restrict the computation of the capture probabilities to $x\in[0:\pi]$. We
perform the following canonical change of variables $(X,x)\rightarrow(z,x)$
with $z=X/G=\cos i$. With these variables, the Hamiltonian (24) becomes

$\mathcal{H}^{\prime}_{\nu_{1}}=-\mathcal{A}Gz^{2}-\mathcal{B}G\left(1-z^{2}\right)\cos
2x+n_{1}z.$ (43)

To compute the libration area, we need to know the equation of the separatrix.
The unstable fixed point (28) is on the separatrix, whose energy is

$\mathcal{H}^{\prime}_{S}=\left(\mathcal{D}-\mathcal{B}\right)G\
,\quad\mathrm{with}\quad\mathcal{D}=\frac{n_{1}^{2}}{4G^{2}\left(\mathcal{A}-\mathcal{B}\right)}\
.$ (44)

The equation of the separatrix is then

$-\mathcal{A}Gz^{2}-\mathcal{B}G\left(1-z^{2}\right)\cos
2x+n_{1}z=\left(\mathcal{D}-\mathcal{B}\right)G\ ,$ (45)

and so, the equations of the superior separatrix, $z_{+}$, and the inferior
one, $z_{-}$, are

$\begin{split}z_{\pm}=&\frac{n_{1}}{2G\left(\mathcal{A}-\mathcal{B}\cos
2x\right)}\\\
&\pm\frac{\sqrt{2\mathcal{B}\sin^{2}x\left(\mathcal{A}-\mathcal{B}\cos
2x-\mathcal{D}\right)}}{\mathcal{A}-\mathcal{B}\cos 2x}\end{split}\ .$ (46)

The libration area encircled by the separatrix is

$\begin{split}A_{lib}&=\int_{0}^{\pi}z_{+}\,dx-\int_{0}^{\pi}z_{-}\,dx\\\
&=2\sqrt{2\mathcal{B}}\int_{0}^{\pi}\frac{\sqrt{\sin^{2}x\left(\mathcal{A}-\mathcal{B}\cos
2x-\mathcal{D}\right)}}{\mathcal{A}-\mathcal{B}\cos 2x}\,dx\\\
&=2\sigma\sqrt{\frac{2}{\tau}}\int_{0}^{\pi}\frac{\sin\left(\frac{x}{2}\right)\sqrt{1-\tau\cos
x}}{1-\sigma\cos x}\,dx\ ,\end{split}$ (47)

with

$\sigma=\frac{\mathcal{B}}{\mathcal{A}}\quad\mathrm{and}\quad\tau=\frac{\sigma}{1-\mathcal{D}\mathcal{A}^{-1}}\
.$ (48)

In order that the separatrix and the libration area exist, we must have
$n_{1}<2G(\mathcal{A}-\mathcal{B})<2G(\mathcal{A}+\mathcal{B})$. Therefore,
$\sigma$ and $\tau$ verify $0<\sigma<\tau<1$. We obtain

$\begin{split}&\int_{0}^{\pi}\frac{\sin\left(\frac{x}{2}\right)\sqrt{1-\tau\cos
x}}{1-\sigma\cos
x}\,dx=\frac{1}{\sigma\sqrt{2}}\left(\sqrt{\tau}\arccos\left(1-\frac{4\tau}{1+\tau}\right)\right.\\\
&\left.-\sqrt{\frac{\tau-\sigma}{1+\sigma}}\arccos\left(1-\frac{4\left(\tau-\sigma\right)}{\left(1-\sigma\right)\left(1+\tau\right)}\right)\right)\
,\end{split}$ (49)

and finally

$\begin{split}A_{lib}=&2\left(\arccos\left(1-\frac{4\tau}{1+\tau}\right)\right.\\\
&\left.-\sqrt{\frac{\tau-\sigma}{\tau\left(1+\sigma\right)}}\arccos\left(1-\frac{4\left(\tau-\sigma\right)}{\left(1-\sigma\right)\left(1+\tau\right)}\right)\right).\end{split}$
(50)

For the total area we get

$\begin{split}A_{tot}&=\int_{0}^{\pi}\,dx-\int_{0}^{\pi}z_{-}\,dx\\\
&=\pi-\int_{0}^{\pi}\frac{n_{1}}{2G\left(\mathcal{A}-\mathcal{B}\cos
2x\right)}dx+\frac{A_{lib}}{2}\\\
&=\pi\left(1-\frac{n_{1}}{2G\sqrt{\mathcal{A}^{2}-\mathcal{B}^{2}}}\right)+\frac{A_{lib}}{2}\
.\end{split}$ (51)

If we consider only the tidal effects on the semi-major axis, in the
expressions of $A_{lib}$ and $A_{tot}$ only the semi-major axis depends on the
time. Then, from expression (41) the capture probability is given by

$P_{\mathrm{cap}}=\frac{\partial A_{lib}}{\partial a}\Big{/}\frac{\partial
A_{tot}}{\partial a}\ ,$ (52)

which does not depend on the dissipation law chosen for the semi-major axis.
From the expressions of $A_{lib}$ (Eq. (50)), and of $A_{tot}$ (Eq. (51)), we
can compute the derivatives with respect to the semi-major axis, and finally
the capture probability.

In Fig. 5, we show the capture probability in the $\nu_{1}$ resonance for
Phobos as a function of the obliquity of Mars (with initial values taken from
Table 1). The solid line gives the theoretical estimation obtained with
expression (52), using a critical semi-major axis obtained with expression
(27). In order to obtain a numerical estimation of the capture probability
with different values of the tidal dissipation (for more details see Eq. (56)
and Sect. 6.2), we perform numerical integrations of the resonant Hamiltonian
(19). We modify this Hamiltonian following Goldreich (1965) and Kinoshita
(1993) to consider the motion of the Martian equator, which we model here with
a constant obliquity and a uniform precession.

For the present tidal dissipation of Mars ($Q=99.5$), there are some
differences between the numerical results and the analytic computation.
However, for weaker tidal dissipation ($Q=5000$), which corresponds to a
slower evolution of the semi-major axis of Phobos (adiabatic evolution), we
observe that there is a good agreement. We hence conclude that the evolution
of Phobos cannot be considered as adiabatic.

Figure 5: Capture probability of Phobos in the $\nu_{1}$ resonance with
respect to the obliquity of Mars. The black curve corresponds to the capture
probability computed with Eq. (52). The blue crosses and the red dots
correspond to the capture probability determined for each value of the
obliquity with $N=720$ numerical integrations of the resonant Hamiltonian (Eq.
(19)), respectively with $Q=99.5$ and $Q=5000$ for the effective specific
tidal dissipation. The errorbars are given by $\pm 1/\sqrt{N}$.

In Fig. 6, we show the theoretical capture probability (Eq. (52)) also as
function of the obliquity of Mars, but for different inclinations and
eccentricities of Phobos. As expected, we observe that the capture probability
increases with the obliquity, because for small obliquity, $\varepsilon$, the
amplitude of the $\nu_{1}$ resonance is proportional to $\varepsilon^{4}$ (Eq.
(19)). On the other hand, we observe that the capture probability strongly
decreases with the inclination of Phobos. The red curve, obtained with $e=0$,
is very close to the black curve, obtained with $e=0.015$, that is, the
capture probability is not much influenced by a small variation in
eccentricity. This was also expected, since the amplitude of the $\nu_{1}$
resonance is proportional to $1+3e^{2}/2$ (Eq. (19)).

Figure 6: Capture probability of Phobos in the $\nu_{1}$ resonance with
respect to the obliquity of Mars for different orbital parameters of Phobos.
The black line corresponds to the current parameters.

The Hamiltonian of the $\nu_{1}$ resonance is similar to the one of the second
order ivection resonance studied by Xu & Fabrycky (2019) in the case of a
binary star system. The resonances called ivection by Xu & Fabrycky (2019)
correspond to a resonance between the node precession frequencies of planets
orbiting around a star with the mean motion of the binary. At the lowest order
in the mutual inclination of planets, Xu & Fabrycky (2019) obtained that the
Hamiltonian of the second order ivection resonance can be written in the same
dimensionless form as the $j+2:j$ orbital resonance described by Borderies &
Goldreich (1984). In the limit of small inclination, it is then also possible
to use the analytical expression for the capture probability given by
Borderies & Goldreich (1984).

The capture probabilities in the $\nu_{2}$ resonance can be obtained with a
similar method as for the $\nu_{1}$ resonance. However, unfortunately it is
not possible to find a simple analytical expression for the corresponding
$A_{lib}$ integral (Eq. (47)). At the lowest order in inclination of the
satellite, the $\nu_{2}$ resonance and the first order ivection resonance of
Xu & Fabrycky (2019) can be written in the same dimensionless form as the
$j+1:j$ orbital resonance described by Borderies & Goldreich (1984). In the
limit of small inclination, the analytical expression of the capture
probability given by Borderies & Goldreich (1984) can then also be used for
the $\nu_{2}$ resonance.

##  5 Interaction between resonances

In previous section, we have studied the dynamics of an isolated resonance.
However, in the full Hamiltonian (17) there is a large number of resonant
terms that may interact. This is particularly true, when two resonant islands
are close to each other, as it is the case of the $\nu_{1}$ resonance and the
evection resonance $\nu_{e}$ (Yokoyama 2002). Indeed, in general we have
$\dot{\varpi}\approx-\dot{\Omega}$, which is the case for a satellite with a
weakly inclined orbit with respect to the equator of its planet. Thus, for
nonzero eccentricities, the libration width of these two resonances may
overlap and trigger chaotic motion (Chirikov 1979). In order to describe this
interaction, we need to add to the previous $\nu_{1}$ resonant Hamiltonian
(19) the term corresponding to the evection resonance $\nu_{e}$. The
interaction Hamiltonian is then

$\begin{split}&\mathcal{H}_{int}=\frac{\mathcal{E}}{2a^{3}\left(1-e^{2}\right)^{3/2}}\left(1-3\cos^{2}i\right)+\frac{3\mathcal{G}m_{0}a^{2}}{2a_{1}^{3}}\left(-\frac{e^{2}}{4}\right.\\\
&+\frac{1}{4}\left(1+\frac{3}{2}e^{2}\right)\left(\sin^{2}i+\sin^{2}\varepsilon-\frac{3}{2}\sin^{2}i\sin^{2}\varepsilon\right)\\\
&-\frac{1}{16}\left(1+\frac{3}{2}e^{2}\right)\sin^{2}i\left(1-\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi+\Omega\right)\right)\\\
&\left.-\frac{5}{64}e^{2}\left(1+\cos
i\right)^{2}\left(1+\cos\varepsilon\right)^{2}\cos\left(2\left(\lambda_{1}-\psi-\varpi\right)\right)\right).\end{split}$
(53)

###  5.1 Frequency analysis method

The interaction Hamiltonian has two degrees of freedom, and it is thus not
integrable. This kind of problem is often studied using surface sections, but
here we propose a different method: we use stability maps based on frequency
analysis (Laskar 1988, 1990; Laskar et al. 1992; Laskar 1993, 2003). This
method decomposes a discrete temporal function in a quasi-periodic
approximation and can estimate its stability. We consider a grid of initial
conditions and integrate the equations of motion over a time $T$. Then, we
perform a frequency analysis222We analyze the quantity
$k^{\prime}_{x}+ik^{\prime}_{y}$, where
$\@vec{k}^{\prime}=(k_{x},k_{y},k_{z})$ corresponds to the coordinates of the
normal to the satellite’s orbit expressed in the referential of the equator of
the planet. In the quasi-periodic approximation, the precession frequency
corresponds to the main frequency with the largest amplitude. over the time
intervals $[0:T/2]$ and $[T/2:T]$, and determine the precession frequency of
the ascending node of the satellite in each interval, $f_{1}$ and $f_{2}$,
respectively. The stability is measured by

$\sigma\equiv\left\lvert 1-\frac{f_{2}}{f_{1}}\right\rvert\ ,$ (54)

which estimates the chaotic diffusion of the precession frequency (Dumas &
Laskar 1993; Laskar 1993). The larger $\sigma$ is, the more unstable the
orbital motion of the satellite is. For stable motion we have $\sigma\sim 0$,
while $\sigma\ll 1$ if the motion is weakly perturbed, and $\sigma\sim 1$ when
the motion is chaotic. It is difficult to know precisely what is the value of
$\sigma$ for which the motion is stable or unstable, but a threshold of
stability $\sigma_{s}$ can be estimated such that most of the trajectories
with $\sigma<\sigma_{s}$ are stable (Couetdic et al. 2010). For
$\sigma<\sigma_{s}$, we still would like to distinguish the circulation
trajectories from those that are in resonance. For that purpose, in those
cases we compute a second quantity

$\theta\equiv\left\lvert 1+\frac{f_{1}}{n^{\prime}_{1}}\right\rvert\ ,$ (55)

which measures the relative difference between the precession frequency and
the resonance frequency $-n^{\prime}_{1}$. We assume that a trajectory is in
resonance when $\theta<\theta_{r}$, where the boundary $\theta_{r}$ is chosen
such that it is able to correctly identify the libration area in the case of
the resonant Hamiltonian (19). In the stability maps, we attribute a color
scheme to the different $\sigma$ going from blue (stable) to red (unstable),
while for the resonant trajectories we apply a black filter.

###  5.2 Stability maps

As for surface sections, the frequency analysis is a numerical method, and
thus we need to attribute values to the parameters of the system. We adopt
here the values for Phobos and Mars (Table 1) near the $\nu_{1}$ resonance,
for which $a_{\nu_{1}}/R\approx 2.617$. This resonance is more interesting
than the $\nu_{2}$ resonance, since it is close to the evection resonance
$\nu_{e}$. We also fix the obliquity of Mars at
$\varepsilon=90\mathrm{\SIUnitSymbolDegree}$, because the resonance width of
$\nu_{1}$ is larger for this obliquity (Eq. (19)) for
$\varepsilon\in[0\mathrm{\SIUnitSymbolDegree}:90\mathrm{\SIUnitSymbolDegree}]$.
We build a 2D mesh of initial conditions where the inclination of the
satellite varies from $0.04\mathrm{\SIUnitSymbolDegree}$ to
$4\mathrm{\SIUnitSymbolDegree}$ with a step size of
$0.04\mathrm{\SIUnitSymbolDegree}$, and the canonical variable
$x=\lambda_{1}-\psi+\Omega$ varies from $0\mathrm{\SIUnitSymbolDegree}$ to
$180\mathrm{\SIUnitSymbolDegree}$ with a step size of
$2\mathrm{\SIUnitSymbolDegree}$, and we integrate the equations of motion over
the time interval $T=80$ kyr.

In Fig. 7, we show the stability maps obtained for different initial values of
the eccentricity of Phobos $e=10^{-5}$, 0.015, 0.030, and 0.045, since the
amplitude of the evection resonance increases with the eccentricity (Eq.
(53)). For each value of the initial eccentricity, we numerically integrate
the equations of motion obtained for the full secular Hamiltonian
$\mathcal{H}_{s}$ (Eq. (5)), for the resonant Hamiltonian
$\mathcal{H}_{\nu_{1}}$ (Eq. (19)), and for the interaction Hamiltonian
$\mathcal{H}_{int}$ (Eq. (53)). For the three cases, we consider that the
Martian obliquity is constant, and that the precession of the Martian equator
is uniform. To consider the motion of the Martian equator, we modify the
resonant Hamiltonian $\mathcal{H}_{\nu_{1}}$ and the interaction Hamiltonian
$\mathcal{H}_{int}$ following Goldreich (1965) and Kinoshita (1993). That is,
we obtain three different maps for each eccentricity value333As the location
of the stable point of the $\nu_{1}$ resonance depends on the eccentricity
(Eq. (27)), and also on the considered Hamiltonian, we need to slightly modify
the value of the semi-major axis in order that the location of the stable
point is located about at the same place for each eccentricity and
Hamiltonian.: the left-hand map corresponds to the integrable resonant
Hamiltonian for the $\nu_{1}$ resonance (Eq. (19)), the middle map corresponds
to the sum of this resonant Hamiltonian with the term due to the evection
resonance (Eq. (53)), and the right-hand map illustrates the true dynamics of
the system obtained with the full secular Hamiltonian (Eq. (5)). The resonant
Hamiltonian (Eq. (19)) gives identical trajectories for the initial conditions
$x$ and $x+\pi$, that is why we only consider $x\in[0^{\circ},180^{\circ}]$.
This symmetry is broken when we introduce the perturbations from the other
resonant terms, but the changes are almost imperceptible.

We attribute a color scheme to the different values of $\sigma$, and estimate
for the stability threshold444Following Couetdic et al. (2010), we estimate
the diffusion $\sigma^{\prime}$ for a shorter interval $[0:T^{\prime}]$, where
$T^{\prime}=T/100$. In general, for a stable trajectory
$\sigma<\sigma^{\prime}$, because the accuracy of the frequency analysis
increases with the size of the interval, while $\sigma>\sigma^{\prime}$ for an
unstable trajectory, as the chaotic diffusion increases with time. We then
construct the histogram of the trajectories verifying $\sigma<\sigma^{\prime}$
with a step of $0.5$. The values of $\sigma$ and $\sigma^{\prime}$ are
averaged over the closest neighbors of the grid of initial conditions to
decrease the effects of chaotic trajectories which have small values of the
diffusion. From this histogram, we determine with a linear interpolation the
value $\log_{10}\sigma_{s}=-4.4$ below which $99\%$ of the trajectories verify
$\sigma<\sigma^{\prime}$., $\log_{10}\sigma_{s}=-4.4$. Therefore, the blue and
green areas correspond to stable circulating trajectories, while the yellow,
orange, and red ones correspond to chaotic motions. In the panels
corresponding to the resonant Hamiltonian (left-hand maps in Fig. 7) we
additionally show the level curves of constant energy obtained with expression
(19), as this problem is integrable (see also Fig. 3). We observe there is a
good agreement between these curves and the stability map obtained with
frequency analysis. These panels thus allow us to calibrate the transition
between the libration and circulation areas. We find that trajectories inside
the separatrix are bounded by $\log_{10}\theta_{r}=-3.9$ (Eq. (55)). We then
apply a black filter for stable trajectories with $\theta<\theta_{r}$ that
correspond to resonant regions555For the maps with $e=0.03$ and $e=0.045$ in
Fig. 7, some black regions appear around $i=0\mathrm{\SIUnitSymbolDegree}$ and
$i=1\mathrm{\SIUnitSymbolDegree}$ that do not correspond to the $\nu_{1}$
resonance. They correspond to stable circulation regions that incidentally
also have $\theta<\theta_{r}$..

resonant Hamiltonian $\mathcal{H}_{\nu_{1}}$ (Eq. (19)) | interaction Hamiltonian $\mathcal{H}_{int}$ (Eq. (53)) |  secular Hamiltonian $\mathcal{H}_{s}$ (Eq. (5))
---|---|---

Figure 7: Stability maps for Phobos (Table 1) using an obliquity for Mars
$\varepsilon=90\mathrm{\SIUnitSymbolDegree}$. From top to bottom, the initial
eccentricity is $e=10^{-5}$, $e=0.015$, $e=0.03$, $e=0.045$. Each column
corresponds to the integration using different Hamiltonians: resonant
Hamiltonian (left), interaction Hamiltonian (middle) and the full secular
Hamiltonian (right). The color scale indicates the value of the quantity
$\log_{10}\sigma$. The blue and green areas correspond to stable circulating
trajectories, while the yellow, orange, and red ones correspond to chaotic
motions. The black areas correspond to stable resonant regions, for which
$\log_{10}\theta<-3.9$.

In Fig. 7, we observe that as we increase the eccentricity, the dynamics
changes considerably between the resonant case and the full problem. However,
a first striking evidence is that the evection resonance $\nu_{e}$ is the
major responsible for this modification, since the middle maps obtained only
with the interaction Hamiltonian (Eq. (53)) are very similar to the right-hand
maps obtained with the full secular Hamiltonian (Eq. (5)). This means that the
dynamics of Phobos around the $\nu_{1}$ resonance is essentially governed by
the interaction between $\nu_{1}$ and $\nu_{e}$. More detailed studies of this
problem, thus only need to take the two-degree of freedom Hamiltonian (53)
into account.

As expected, for a nearly zero eccentricity, $e=10^{-5}$ (first row in Fig.
7), all maps are similar to the one obtained with the resonant Hamiltonian
(19). Indeed, for such a small eccentricity, the amplitude of the evection
resonance $\nu_{e}$ is nearly zero, and the dynamics is dominated by the
$\nu_{1}$ resonance. Nonetheless, in the case of the full Hamiltonian we
observe already some diffusion of the orbits in the region around the
libration area, because the eccentricity is not exactly zero. For the present
eccentricity of Phobos, $e=0.015$ (second row in Fig. 7), the $\nu_{1}$
resonance still presents a large libration area in the full Hamiltonian case,
but is now surrounded by a large chaotic region. Its chaotic nature is notably
characterized by a high fickleness of the value of the diffusion index,
$\sigma$. For an eccentricity of $e=0.03$ (third row in Fig. 7), the libration
area corresponding to the $\nu_{1}$ resonance still exists, but it is
considerably smaller and deformed, while for $e=0.045$ (last row in Fig. 7)
the $\nu_{1}$ resonance almost disappeared in the middle of a chaotic region.
Finally, we also observe for the full secular Hamiltonian and the resonant
Hamiltonian that as the eccentricity increases, the center of the libration
area is slightly shifted toward increasing values of $x>90^{\circ}$.

In Fig. 7, we confirm that the interaction between the resonances $\nu_{1}$
and $\nu_{e}$ increases with the eccentricity. For the present eccentricity of
Phobos (or smaller), this interaction is not very important. As a result, in
the future Phobos may survive in the $\nu_{1}$ resonance. However, as the
eccentricity increases, we observe that the chaotic region around $\nu_{1}$
progressively replaces the libration area. It is thus unlikely for a satellite
with eccentricities larger than $0.05$ to be captured in this resonance.

##  6 Application to Phobos

In order to maximize the impact of the eviction-like resonances on the
inclination of a satellite, the semi-major axis must decrease. Tidal
interactions between the satellite and the planet can account for this
evolution, provided that the orbital period of the satellite is shorter than
the rotational period of the planet, which is the case of Phobos, the largest
satellite of Mars (e.g., Szeto 1983). Among the main satellites of solar
system planets, only Triton, the largest satellite of Neptune, is also
spiraling into the planet, but because it is on a retrograde orbit (e.g.,
Correia 2009). Therefore, in this section we apply our model to study Phobos.

At present, the semi-major axis of Phobos is $a/R=2.761$, and the eccentricity
is $e=0.01511$ (Jacobson & Lainey 2014). Among the fourteen resonances listed
in Eq. (18), seven resonances can act on a prograde satellite such as Phobos.
During its past evolution, Phobos already encountered four of these
resonances, but their effect on its orbit is presumably to have been
negligible (Yokoyama 2002), since the amplitudes of these resonances are
proportional to the square of the eccentricity. In the future, Phobos will
encounter the $\nu_{e}$ and the $\nu_{1}$ resonances, which occur almost
simultaneously around $a_{\nu_{1}}/R=2.617$, and later the $\nu_{2}$
resonance, which is located at $a_{\nu_{2}}/R=2.147$. For the $\nu_{1}$ and
the $\nu_{2}$ resonances, the amplitudes are not null for eccentricities close
to zero, and so they can modify the future inclination of Phobos.

In Sect. 4.3, we have obtained the capture probability with an analytical
computation. However, expression (52) is only valid for the $\nu_{1}$
resonance and if the evolution is adiabatic. Moreover, it does not take the
interaction between the $\nu_{1}$ resonance and the evection resonance
$\nu_{e}$ into account, that we have seen in Sect. 5. It is then necessary to
perform numerical simulations to capture the correct dynamics while crossing
the $\nu_{1}$ and $\nu_{2}$ resonances. Here we investigate the capture
probabilities in these two resonances, and determine which parameters can
influence this capture. In the case where capture is not possible, we also
study the effects due to the crossing of these resonances on the orbital
evolution of Phobos.

###  6.1 Numerical model

We numerically integrate the set of equations (7) to obtain the secular
orbital evolution of Phobos forced by the orbital motion of Mars and its
equator. We consider that the orbit of Mars in the ecliptic is circular and
that its equator has a uniform precession and a constant obliquity.

Tides raised by Phobos on Mars, modify all the orbital parameters of the
satellite with time. However, for simplicity we consider only the evolution of
the semi-major axis. We also do not consider tides raised by Mars on Phobos,
whose effect is very small, since the rotation of Phobos is currently
synchronous and on a near circular orbit. In order to describe the evolution
of the semi-major axis, we adopt a constant$-Q$ tidal model (e.g., Szeto
1983),

$\dot{a}=-3\sqrt{\frac{\mathcal{G}}{m_{1}}}\frac{mR^{5}}{a^{11/2}}\frac{k_{2}}{Q},$
(56)

where $m$ is the mass of Phobos, $R$ is the radius of Mars, $k_{2}$ is the
second Love number, and $Q$ is the effective specific tidal dissipation. From
the observation of the satellites of Mars, it is possible to put constraints
on the coefficients $k_{2}$ and $Q$. We adopt the values obtained by Jacobson
& Lainey (2014) with this method: $k_{2}=0.183$ and $Q=99.5$ (Table 1).

###  6.2 Numerical setup

The capture probability in eviction-like resonances depends on the inclination
of the satellite and on the obliquity of the planet (Sect. 4.3). However, it
can also depend on the eccentricity, because of the interaction with neighbor
resonances (Sect. 5). The current eccentricity of Phobos is about $0.015$, but
it is still being damped to smaller values. It is then likely that it will be
smaller when Phobos approaches the resonances. Then, in order to observe the
impact of the eccentricity in the capture probability, we adopt different
eccentricity values up to 0.015.

The obliquity of Mars is chaotic and cannot be precisely computed beyond about
10 Myr (Touma & Wisdom 1993; Laskar & Robutel 1993; Laskar et al. 2004). It
experiences variations between $0\mathrm{\SIUnitSymbolDegree}$ and
$60\mathrm{\SIUnitSymbolDegree}$ over 50 Myr, and can even reach values beyond
$60\mathrm{\SIUnitSymbolDegree}$ over higher periods of time. For instance,
over 100 Myr, the obliquity has a probability of $7.23\%$ to reach
$60\mathrm{\SIUnitSymbolDegree}$ and cannot have values beyond
$70\mathrm{\SIUnitSymbolDegree}$, but over 5 Gyr, the probabilities for the
obliquity to reach $60\mathrm{\SIUnitSymbolDegree}$,
$70\mathrm{\SIUnitSymbolDegree}$, and $80\mathrm{\SIUnitSymbolDegree}$, are
respectively $95.35\%$, $8.51\%$, and $0.015\%$ (Laskar et al. 2004). Due to
the decrease in the semi-major axis, it is expected that Phobos will collide
with Mars in about 40 Myr (Efroimsky & Lainey 2007). It is then very unlikely
that Mars has an obliquity higher than $60\mathrm{\SIUnitSymbolDegree}$ when
Phobos encounters the eviction-like resonances. We nevertheless investigate
the probability of capture for all obliquity values up to
$90\mathrm{\SIUnitSymbolDegree}$, in order to get a global vision of this
mechanism.

To numerically estimate the capture probability in resonance, for each initial
value of the eccentricity and obliquity we perform $720$ simulations varying
the initial longitude of the ascending node $\Omega$ from
$0\mathrm{\SIUnitSymbolDegree}$ to $359.5\mathrm{\SIUnitSymbolDegree}$ with a
step of $0.5\mathrm{\SIUnitSymbolDegree}$, while fixing the remaining initial
parameters.

The rotational and orbital motion of Mars is complex, but can be approximated
with quasi-periodic solutions (e.g., Ward 1979; Laskar 1990; Laskar & Robutel
1993; Laskar et al. 2004). In such solutions, the Martian obliquity is not
constant, but the current dominant term corresponds to rapid oscillations,
which can induce variations between $10.8\mathrm{\SIUnitSymbolDegree}$ and
$38.0\mathrm{\SIUnitSymbolDegree}$ with a period of about 120 kyr (Ward 1979).
Therefore, its evolution time-scale is much longer than the precession motion
of Phobos, and so it can be considered constant during the resonance crossing.
We consider a fixed orbit for Mars and a uniform precession of the Martian
equator, that is, without taking the planetary perturbations into account that
induce secular variations in periods much longer than the precession motion of
Phobos. We also consider a circular orbit for Mars to simplify the problem.

###  6.3 Resonance $\nu_{1}$

####  6.3.1 Capture probabilities

We compute the capture probabilities in the $\nu_{1}$ resonance using the
present inclination of Phobos, $i=1\mathrm{\SIUnitSymbolDegree}$. For the
eccentricity, we adopt values $e=0$, $10^{-5}$, $10^{-4}$, $10^{-3}$ and
$0.015$, while for the obliquity we adopt values from
$0\mathrm{\SIUnitSymbolDegree}$ to $60\mathrm{\SIUnitSymbolDegree}$ with a
step of $10\mathrm{\SIUnitSymbolDegree}$, and from
$60\mathrm{\SIUnitSymbolDegree}$ to $90\mathrm{\SIUnitSymbolDegree}$ with a
step of $2\mathrm{\SIUnitSymbolDegree}$, since the variations in the capture
probability are more important when the obliquity is high (Sect. 4.3). We
start our simulations for a semi-major axis slightly larger than $a_{\nu_{1}}$
(Eq. (30)) and integrate the equations of motion for 500 kyr. The results are
shown in Fig. 8.

As expected, we observe that the capture probability in the $\nu_{1}$
resonance increases with the obliquity (Sect. 4.3). For all the considered
eccentricities, we did not observe any capture for an obliquity smaller than
$60\mathrm{\SIUnitSymbolDegree}$, and when
$\varepsilon=60\mathrm{\SIUnitSymbolDegree}$, a unique capture occurred (for
$e=10^{-4}$). The probability of capture for $e=0$ increases from $0.28\%$ for
$\varepsilon=66\mathrm{\SIUnitSymbolDegree}$ to $81.4\%$ for
$\varepsilon=90\mathrm{\SIUnitSymbolDegree}$. For an initial eccentricity of
$0.015$, which corresponds to the current value of the eccentricity, the
probability of capture is zero except for
$\varepsilon=82\mathrm{\SIUnitSymbolDegree}$ and
$\varepsilon=90\mathrm{\SIUnitSymbolDegree}$, where one and two captures have
been observed, which correspond to capture probabilities of about $0.14\%$ and
$0.28\%$, respectively. These results are consistent with those observed by
Yokoyama (2002), but present some disagreement with the theoretical
estimations done in Sect. 4.3. One reason is because tides acting on Phobos
semi-major axis are very strong (Eq. (56)), and hence the evolution is not
adiabatic (see Fig. 5). The other reason is the perturbations from the
evection resonance (Sect. 5).

Figure 8: Capture probabilities of Phobos in the $\nu_{1}$ resonance with
respect to the obliquity of the equator of Mars for different initial
eccentricities.

According to expression (19), the influence of the eccentricity should be
weak, because the resonance width slightly increases with the eccentricity.
However, in Fig. 8, we note that the probability of capture sharply decreases
when the eccentricity slightly increases. For instance, for an obliquity of
$90\mathrm{\SIUnitSymbolDegree}$, the probability is $81.4\%$ for $e=0$,
$61.7\%$ for $e=10^{-5}$, $19.2\%$ for $e=10^{-4}$, and $3.6\%$ for
$e=10^{-3}$. This is somehow surprising, because in Fig. 7 we observe that for
the present eccentricity ($e=0.015$) a large resonance island is still
present. However, as seen in previous sections, the $\nu_{1}$ resonance is
very close to the evection resonance $\nu_{e}$. As the semi-major axis of
Phobos is decreasing, this resonance occurs almost simultaneously with the
$\nu_{1}$ resonance. Capture in the evection resonance is not possible, but
crossing it induces a variation in the eccentricity of Phobos as we can see in
Figs. 9, 10, and 11.

Figure 9: Evolution of the inclination (a) and of the eccentricity (b) during
the capture in the $\nu_{1}$ resonance for an obliquity of
$\varepsilon=90\mathrm{\SIUnitSymbolDegree}$ with the initial conditions
$e=0.015$ and $\Omega=351.5\mathrm{\SIUnitSymbolDegree}$. The blue curve
corresponds to the stable fixed point of the $\nu_{1}$ resonance computed with
Eq. (27).  Figure 10: Evolution of the inclination (a) and of the eccentricity
(b) during the crossing of the $\nu_{1}$ resonance for an obliquity of
$\varepsilon=60\mathrm{\SIUnitSymbolDegree}$ with the initial conditions
$e=0.015$, and $\Omega=231\mathrm{\SIUnitSymbolDegree}$.  Figure 11: Evolution
of the inclination (a) and of the eccentricity (b) during the crossing of the
$\nu_{1}$ resonance for an obliquity of
$\varepsilon=60\mathrm{\SIUnitSymbolDegree}$ with the initial conditions
$e=0.015$, and $\Omega=203\mathrm{\SIUnitSymbolDegree}$.

In Fig. 12, we plot the maximum and minimum differences, $\Delta e$, between
the mean eccentricity before and after the encounter with the $\nu_{1}$
resonance as a function of the obliquity, when the capture succeeds and fails.
If the capture fails, regardless of the initial value of the eccentricity
before encountering the $\nu_{1}$ resonance, the final mean eccentricity
always increases. We note that during the crossing of the $\nu_{1}$ resonance
the eccentricity can even take higher values before decreasing to a smaller
final value (Figs. 10, 11). For small obliquities, the eccentricity is always
strongly excited by about $0.04$ (e.g., for
$\varepsilon=0\mathrm{\SIUnitSymbolDegree}$, we have $0.0365<\Delta
e<0.0395$), which is consistent with the results previously observed by Yoder
(1982) and Yokoyama (2002). For obliquities larger than
$60\mathrm{\SIUnitSymbolDegree}$, we observe a smaller average variation
around $0.02$, that can nevertheless still temporarily increase to about
$0.04$ (e.g., for $\varepsilon=90\mathrm{\SIUnitSymbolDegree}$, we have
$0.005<\Delta e<0.0298$). For these larger eccentricity values, the libration
width of the $\nu_{1}$ resonance is much smaller and the capture probability
is hence considerably reduced (see Fig. 7).

In Fig. 12, we also show the statistics when capture succeeds666We do not have
many examples for the initial eccentricity $0.015$ after capture, because the
probability is very small.. In this case, we observe smaller variations in the
eccentricity. Therefore, the chances of capture in the $\nu_{1}$ resonance are
weak because of the nearby evection resonance $\nu_{e}$. This resonance leads
to an increase in the eccentricity, which in turn gives rise to a chaotic
region around the $\nu_{1}$ resonance that reduces its libration width (see
Sect. 5).

Figure 12: Evolution of the maximum and minimum variations in the mean
eccentricity $\Delta e$ owing to the encounter with the $\nu_{1}$ resonance in
the cases where the capture in the resonance fails (a), and succeeds (b), with
respect to the obliquity for different initial eccentricities.

We conclude that capture in the $\nu_{1}$ resonance can only occur for very
small eccentricity values and for an obliquity higher than
$60\mathrm{\SIUnitSymbolDegree}$. These chances are on the low side, because
as we have just seen, the eccentricity of Phobos will likely increase after
crossing the evection resonance, and the obliquity of Mars is not expected to
exceed $60\mathrm{\SIUnitSymbolDegree}$ in the next 40 Myr. It is then
extremely unlikely that in the future Phobos gets caught in the $\nu_{1}$
resonance. Nevertheless, if ever capture occurs, this resonance is stable and
the inclination can increase. Indeed, in Fig. 9 we show an example of this
capture and subsequent evolution for $e=0.015$ and $\varepsilon=90^{\circ}$.
We also confirm that, when capture occurs, the inclination follows the
equilibrium point given by expression (27).

####  6.3.2 Inclination variation

When capture in resonance does not occur, we still observe inclination
variations during the resonance crossing. As for the capture probability, the
exact behavior of the inclination depends on the initial value of the
longitude of the ascending node. In Figs. 10 and 11, we show two examples of
such variations in the case of the $\nu_{1}$ resonance, which correspond to an
increase and to a decrease in the initial inclination, respectively. These
variations modify the state of a satellite and can lead to misinterpretations
of its origin. We can profit from the simulations already performed to
estimate the capture probabilities to put limits on the inclination
variations. Moreover, the $\nu_{2}$ resonance is encountered after the
$\nu_{1}$ resonance, and so it is important to estimate the inclination of
Phobos after crossing the $\nu_{1}$ resonance, because the capture probability
depends on the inclination (Sect. 4.3).

In Fig. 13 we show the maximum and minimum differences between the mean
inclination of Phobos before and after the crossing of the $\nu_{1}$
resonance. As for the capture probability, we observe that these variations
depend on the initial eccentricity and obliquity. The variations increase in
general with the obliquity: they are nearly zero for obliquities below
$40\mathrm{\SIUnitSymbolDegree}$, but they can present an amplitude variation
of about $1\mathrm{\SIUnitSymbolDegree}$ for obliquities larger than
$60^{\mathrm{\SIUnitSymbolDegree}}$. This can seem a modest variation, but
given that the present inclination of Phobos is also
$1\mathrm{\SIUnitSymbolDegree}$, it corresponds to a change around $100\%$. If
the eccentricity of Mars is considered, Yokoyama et al. (2005) observe that,
for some initial conditions, the eccentricity and the inclination after the
crossing of the resonance can even increase to about $0.08$ and
$4\mathrm{\SIUnitSymbolDegree}$, respectively.

Figure 13: Evolution of the maximum and minimum variations in the mean
inclination $\Delta i$ owing to the encounter with the $\nu_{1}$ resonance in
the case where the capture in the resonance fails, with respect to the
obliquity for different initial eccentricities.

###  6.4 Resonance $\nu_{2}$

####  6.4.1 Capture probabilities

As the obliquity of Mars is not expected to exceed
$60\mathrm{\SIUnitSymbolDegree}$ in the next 40 Myr, capture in the $\nu_{1}$
resonance is very unlikely, and so the inclination of Phobos is expected to be
still close to $1\mathrm{\SIUnitSymbolDegree}$ at the moment of the encounter
with the $\nu_{2}$ resonance (Sect. 6.3.2).

We compute the capture probabilities in the $\nu_{2}$ resonance using the
present inclination of Phobos, $i=1\mathrm{\SIUnitSymbolDegree}$. For the
eccentricity, we adopt two values, $e=0$ and $0.015$, while for the obliquity
we adopt values varying from $0\mathrm{\SIUnitSymbolDegree}$ to
$90\mathrm{\SIUnitSymbolDegree}$ with a step of
$10\mathrm{\SIUnitSymbolDegree}$. We start our simulations for a semi-major
axis slightly larger than $a_{\nu_{2}}$ (Eq. (40)) and integrate the equations
of motion for 250 kyr. The capture probability statistics are shown in Fig.
14.

Figure 14: Probabilities of capture of Phobos in the $\nu_{2}$ resonance with
respect to the obliquity of the equator of Mars for different initial
eccentricities and inclinations.

We observe that, as for the $\nu_{1}$ resonance, the probability of capture
increases with the obliquity. It is zero for an obliquity smaller than
$10\mathrm{\SIUnitSymbolDegree}$ and $100\%$ if the obliquity is higher than
$30\mathrm{\SIUnitSymbolDegree}$. These results are also consistent with those
observed by Yokoyama (2002). For an obliquity of
$20\mathrm{\SIUnitSymbolDegree}$, the probability is larger for a null initial
eccentricity than for the current eccentricity, although the effect of the
eccentricity is not very important. Indeed, contrary to the $\nu_{1}$
resonance, which is perturbed by the nearby evection resonance, the $\nu_{2}$
resonance is isolated from the remaining resonances, and so the eccentricity
only slightly modifies the amplitude of the libration width (Eq.
(LABEL:Ham_nu2)).

In Fig. 15, we show an example for the evolution of the eccentricity and
inclination of Phobos in the case where the capture succeeds. We confirm that
the inclination follows the equilibrium point given by expression (37) and
that the eccentricity remains almost unchanged.

Figure 15: Evolution of the inclination (a), and of the eccentricity (b)
during the capture in the $\nu_{2}$ resonance for an obliquity of
$\varepsilon=20\mathrm{\SIUnitSymbolDegree}$ with the initial conditions
$e=0.015$ and $\Omega=0\mathrm{\SIUnitSymbolDegree}$. The blue curve
corresponds to the stable fixed point of the $\nu_{2}$ resonance computed with
Eq. (37).

The present obliquity of Mars is about $25\mathrm{\SIUnitSymbolDegree}$ and
should have variations between $10\mathrm{\SIUnitSymbolDegree}$ and
$40\mathrm{\SIUnitSymbolDegree}$ for the next 10 Myr (Laskar et al. 2004).
Contrary to the $\nu_{1}$ resonance, we thus conclude that in the future
Phobos is likely to be captured in the $\nu_{2}$ resonance. As a result, we
expect that the inclination of Phobos increases during the last stages of its
evolution following the equilibrium point given by expression (39). Therefore,
although Phobos spends most of its life near the equatorial plane of Mars, we
cannot rule out that the impact between Phobos and the surface of Mars will
occur at high latitudes (see Fig. 4). Yokoyama et al. (2005) note, however,
that escape from the $\nu_{2}$ resonance is possible from an inclination of
$33\mathrm{\SIUnitSymbolDegree}$ due to interactions with several resonances.

The modification observed in the inclination during the crossing of the
$\nu_{1}$ resonance without capture (Sect. 6.3.2) can impact the capture
probability in the subsequent $\nu_{2}$ resonance. As the capture probability
decreases when the inclination increases (Sect. 4.3), Yokoyama et al. (2005)
noted that for $i=4\mathrm{\SIUnitSymbolDegree}$ capture in the $\nu_{2}$
resonance is uncertain. Therefore, in Fig. 14 we also compute the capture
probability for $i=4\mathrm{\SIUnitSymbolDegree}$ (blue curve) as a function
of the obliquity. Indeed, we observe that the capture probability decreases
for obliquities smaller than $60\mathrm{\SIUnitSymbolDegree}$, but for
$20\mathrm{\SIUnitSymbolDegree}<\varepsilon<40\mathrm{\SIUnitSymbolDegree}$
the capture is still possible, although with a smaller probability around
$20\%$. However, for obliquities larger than $70\mathrm{\SIUnitSymbolDegree}$
the capture is certain even for $i=4\mathrm{\SIUnitSymbolDegree}$.

####  6.4.2 Inclination variation

In Fig. 16, we show the evolution of the mean inclination during the crossing
of the $\nu_{2}$ resonance when the capture does not occur. We only plot the
results for the obliquities $0\mathrm{\SIUnitSymbolDegree}$,
$10\mathrm{\SIUnitSymbolDegree}$, and $20\mathrm{\SIUnitSymbolDegree}$,
because for values of the obliquity larger than
$30\mathrm{\SIUnitSymbolDegree}$ the capture always occurs. As for the
$\nu_{1}$ resonance, the amplitude of the variation increases with the
obliquity. However, we did not observe any significant variation in the
eccentricity during the crossing of the $\nu_{2}$ resonance.

Figure 16: Evolution of the maximum and minimum variations $\Delta i$ in the
mean inclination owing to the encounter with the $\nu_{2}$ resonance in the
case where the capture fails, with respect to the obliquity for different
initial eccentricities.

##  7 Conclusion

In this paper, we have revisited the secular dynamics of a satellite while
perturbed by the orbital motion of its planet. Previous studies have shown
that some of these resonances can significantly modify the eccentricity
(evection-like resonance) or the inclination (eviction-like resonance) of a
satellite (e.g., Yoder 1982; Touma & Wisdom 1998; Yokoyama 2002; Yokoyama et
al. 2005). Eviction-like resonances are important to map, because the
inclination can be used to put constrains on the formation of the satellite.
Moreover, some eviction-like resonances can be present even for the nearly
zero eccentricities currently observed for the main satellites of the solar
system planets.

Capture in eviction-like resonances is a possibility when the semi-major axis
of the planet is decreasing. This capture can lead to a significant increase
in the inclination of the satellite. However, for nearly zero eccentricity,
capture is only possible in two of these resonances that we dubbed $\nu_{1}$
and $\nu_{2}$. These resonances are placed close to the planets’ surface, and
only satellites spiraling into the planet may be affected, such as Phobos or
Triton. We analytically compute the capture probability in the $\nu_{1}$
resonance and compare it to numerical simulations. We have shown that the
evolution of Phobos cannot be considered as adiabatic.

When the eccentricity of the satellite is nonzero, evection-like and
additional eviction-like resonances are possible. Some of these resonances are
close to each other, and they can interact between them. We have studied these
interactions with the method of frequency map analysis that can be used as an
alternative to surface sections, since it is easier to implement in problems
with many degrees of freedom. In the case of the $\nu_{1}$ resonance, we
confirm that the interaction with the evection resonance $\nu_{e}$ introduces
a chaotic region around the $\nu_{1}$ resonance (Yokoyama 2002; Yokoyama et
al. 2005). In the case of Phobos, for eccentricities larger than 0.05, the
evection resonance completely destroys the stability of the $\nu_{1}$
resonance and capture is no longer possible. In addition, because the evection
resonance occurs almost simultaneously with the $\nu_{1}$ resonance, a
satellite that is spiraling into the planet will have its eccentricity excited
by the evection resonance, which in turn decreases the chances of capture in
the $\nu_{1}$ eviction-like resonance.

When the orbit of the satellite crosses a given secular resonance, capture is
not certain. For eviction-like resonances, the probability of capture
increases with the obliquity of the planet and decreases with the inclination
of the satellite. When capture fails, the inclination of the satellite still
undergoes some changes, although much smaller than in the capture scenario,
and it can either increase or decrease. In the case of the $\nu_{1}$
resonance, the capture probability also strongly depends on the eccentricity.
However, this is a collateral effect that results from the proximity of the
nearby evection resonance.

We apply our model to the study of Phobos, whose orbit will be modified in the
future due to the encounters with the $\nu_{1}$ and $\nu_{2}$ eviction-like
resonances. We obtain for the capture probabilities similar results to
Yokoyama (2002). The capture in the $\nu_{1}$ resonance is almost impossible,
given the present obliquity of Mars and the present eccentricity of Phobos.
Nevertheless, the crossing of this resonance can lead to variations in the
inclination between $-0.08\mathrm{\SIUnitSymbolDegree}$ and
$0.07\mathrm{\SIUnitSymbolDegree}$. On the contrary, capture in the following
$\nu_{2}$ resonance is almost certain if the obliquity of Mars and the
inclination of Phobos are similar to the current values, which will lead to a
continuous increase in the inclination. If the capture fails, the crossing of
the $\nu_{2}$ resonance can lead to variations in the inclination between
$-0.7\mathrm{\SIUnitSymbolDegree}$ and $0.4\mathrm{\SIUnitSymbolDegree}$.

We have seen that secular resonances due to stellar perturbations, such as
evection-like and eviction-like resonances, can significantly influence the
motion of a satellite. The effect of planetary perturbations on satellite
orbits are often considered as negligible, but they also induce secular
resonances, whose crossings can lead to a one-time variation in the
eccentricity or the inclination. If the perturbations are strong enough,
captures in these resonances could also be possible, and induce more
significant changes. Future studies on the orbital evolution of satellite
orbits should thus also estimate the effect from these additional secular
resonances.

###### Acknowledgements.

This work was supported by CIDMA (UIDB/04106/2020 and UIDP/04106/2020), CFisUC
(UIDB/04564/2020 and UIDP/04564/2020), PHOBOS (POCI-01-0145-FEDER-029932),
ENGAGE SKA (POCI-01-0145-FEDER-022217), and GRAVITY (PTDC/FIS-AST/7002/2020),
funded by COMPETE 2020 and FCT, Portugal.

## References

  * Agnor & Hamilton (2006) Agnor, C. B. & Hamilton, D. P. 2006, Nature, 441, 192
  * Batygin & Morbidelli (2020) Batygin, K. & Morbidelli, A. 2020, ApJ, 894, 143
  * Borderies & Goldreich (1984) Borderies, N. & Goldreich, P. 1984, Celestial Mechanics, 32, 127
  * Boué & Laskar (2006) Boué, G. & Laskar, J. 2006, Icarus, 185, 312
  * Canup (2005) Canup, R. M. 2005, Science, 307, 546
  * Canup & Asphaug (2001) Canup, R. M. & Asphaug, E. 2001, Nature, 412, 708
  * Carvalho & Vilhena de Moraes (2020) Carvalho, J. P. S. & Vilhena de Moraes, R. 2020, European Physical Journal Special Topics, 229, 1491
  * Chandrasekhar (1987) Chandrasekhar, S. 1987, Ellipsoidal figures of equilibrium
  * Chirikov (1979) Chirikov, B. V. 1979, Physics Reports, 52, 263
  * Correia (2009) Correia, A. C. M. 2009, ApJ, 704, L1
  * Couetdic et al. (2010) Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10
  * Dumas & Laskar (1993) Dumas, H. S. & Laskar, J. 1993, Phys. Rev. Lett., 70, 2975
  * Efroimsky & Lainey (2007) Efroimsky, M. & Lainey, V. 2007, Journal of Geophysical Research (Planets), 112, 12003
  * Farago & Laskar (2010) Farago, F. & Laskar, J. 2010, MNRAS, 401, 1189
  * Frouard et al. (2010) Frouard, J., Fouchard, M., & Vienne, A. 2010, A&A, 515, A54
  * Goldreich (1965) Goldreich, P. 1965, AJ, 70, 5
  * Henrard (1982) Henrard, J. 1982, Celestial Mechanics, 27, 3
  * Henrard (1993) Henrard, J. 1993, The Adiabatic Invariant in Classical Mechanics, ed. C. K. R. T. Jones, U. Kirchgraber, & H. O. Walther (Berlin, Heidelberg: Springer Berlin Heidelberg), 117–235
  * Inderbitzi et al. (2020) Inderbitzi, C., Szulágyi, J., Cilibrasi, M., & Mayer, L. 2020, MNRAS, 499, 1023
  * Jacobson (2010) Jacobson, R. A. 2010, AJ, 139, 668
  * Jacobson & Lainey (2014) Jacobson, R. A. & Lainey, V. 2014, Planet. Space Sci., 102, 35
  * Kinoshita (1993) Kinoshita, H. 1993, Celestial Mechanics and Dynamical Astronomy, 57, 359
  * Konopliv et al. (2016) Konopliv, A. S., Park, R. S., & Folkner, W. M. 2016, Icarus, 274, 253
  * Lambeck (1979) Lambeck, K. 1979, J. Geophys. Res., 84, 5651
  * Laskar (1988) Laskar, J. 1988, A&A, 198, 341
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Laskar (1993) Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257
  * Laskar (2003) Laskar, J. 2003, arXiv Mathematics e-prints, math/0305364
  * Laskar et al. (2004) Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004, Icarus, 170, 343
  * Laskar et al. (1992) Laskar, J., Froeschlé, C., & Celletti, A. 1992, Physica D Nonlinear Phenomena, 56, 253
  * Laskar & Robutel (1993) Laskar, J. & Robutel, P. 1993, Nature, 361, 608
  * Nesvorný et al. (2007) Nesvorný, D., Vokrouhlický, D., & Morbidelli, A. 2007, AJ, 133, 1962
  * Smart (1953) Smart, W. M. 1953, Celestial Mechanics (London, New York, Longmans, Green)
  * Szeto (1983) Szeto, A. M. K. 1983, Icarus, 55, 133
  * Touma & Wisdom (1993) Touma, J. & Wisdom, J. 1993, Science, 259, 1294
  * Touma & Wisdom (1998) Touma, J. & Wisdom, J. 1998, AJ, 115, 1653
  * Tremaine et al. (2009) Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706
  * Ward (1979) Ward, W. R. 1979, J. Geophys. Res., 84, 237
  * Xu & Fabrycky (2019) Xu, W. & Fabrycky, D. 2019, arXiv e-prints, arXiv:1904.02290
  * Yoder (1979) Yoder, C. F. 1979, Celestial Mechanics, 19, 3
  * Yoder (1982) Yoder, C. F. 1982, Icarus, 49, 327
  * Yokoyama (2002) Yokoyama, T. 2002, Planet. Space Sci., 50, 63
  * Yokoyama et al. (2005) Yokoyama, T., Mana, M. R., do Nascimento, C., Santos, M. T., & Callegari, N., J. 2005, A&A, 429, 731