## 2 Theoretical Background
### The Disturbing Functions  
Following [12], the solar disturbing function can be written as:  
\[ \begin{split}&\mathcal{R}_{\odot}=\sum_{l=2}^{\infty}\sum_{m,p,q=0} ^{l}\sum_{j,r=-\infty}^{\infty}\mu_{\odot}\big{(}\frac{a^{l}}{d_{\odot}^{l+1}} \big{)}e_{m}\frac{(l-m)^{3}}{(l+m)^{3}}F_{lmp}(i)F_{lmq}(i_{\odot})H_{lpy}(e)G_ {lqr}(e_{\odot})\times\\ &\cos\big{[}(l-2p+j)M-(l-2q+r)M_{\odot}+(l-2p)\omega-(l-2q)\omega _{\odot}+m(\Omega-\Omega_{\odot})\big{]},\end{split}  \tag{1}\]  
where \(\mu_{\odot}\) is the gravitational parameter of the Sun, the Keplerian elements of both the satellite and the Sun are written with respect to the equatorial reference plane, \(F_{lmp}\) and \(F_{lmq}\) are [PERSON]'s inclination functions, \(H_{lpy}\) and \(G_{lqr}\) are [PERSON] coefficients and \(\epsilon_{m}\) can be found in [12]. On the other hand, the motion of the Moon around the Earth is significantly perturbed by the gravitational attraction of the Sun; hence, the lunar inclination, node and argument of the perigee with respect to the celestial equator evolve as nonlinear functions of time. To overcome this issue, usually, a mixed reference plane is adopted, where the elements of the satellite are written with respect to the equatorial plane, while the elements of the Moon are referred to as the ecliptic, so that \(i_{\mathcal{Q}}\) is approximately constant and \(\omega_{\mathcal{Q}}\) and \(\Omega_{\mathcal{Q}}\) are approximately linear functions of time [12]. In this way, the lunar disturbing function reads  
\[ \begin{split}&\mathcal{R}_{\mathcal{Q}}\ =\sum_{l=2}^{\infty}\sum_{m,p,s,q=0}^{l +\infty}\sum_{j,r=-\infty}^{+\infty}(-1)^{l+s}(-1)^{l+}\frac{\mu_{\mathcal{Q} }^{l}}{2a_{\mathcal{Q}}}\frac{(l-s)!}{(l+m)!}(\frac{a}{a_{\mathcal{Q}}})^{l} F_{lmp}(i)F_{lsq}(i_{\mathcal{Q}})\times\\ & H_{lpy}(e)G_{lqr}(e_{\mathcal{Q}}\ \big{\{}(-1)^{k_{2}}U_{l}^{m,-s}(e) \cos[(l-2p+j)M+(l-2q+r)M_{\mathcal{Q}}\ +\\ &(l-2p)\omega+(l-2q)\omega_{\mathcal{Q}}\ +m\Omega+s(\Omega_{\mathcal{Q}}\ -\frac{\pi}{2})-y_{s}\pi]\\ &+(-1)^{k_{3}}U_{l}^{m,s}(\epsilon)\cos[(l-2p+j)M-(l-2q+r)M_{ \mathcal{Q}}\ +\\ &(l-2p)\omega-(l-2q)\omega_{\mathcal{Q}}\ +m\Omega-s(\Omega_{\mathcal{Q}}\ -\frac{\pi}{2})-y_{s}\pi]\Big{\}},\end{split}  \tag{2}\]  
where \(\mu_{\mathcal{Q}}\) is the gravitational parameter of the Moon, \(\epsilon\) is the obliquity of the ecliptic, \(y_{s}=0\) for \(s\) even while \(y_{s}=\frac{1}{2}\) for \(s\) odd. The analytical expansions of the Kaula inclination functions and of \(U_{m}^{l,\pm s}\), the Hansen coefficients and the coefficients \(\epsilon_{m}\), \(\epsilon_{s}\), \(k_{1}\), \(k_{2}\), \(k_{3}\) can be found in [12; 13; 14]. We are interested in a model including the oblateness effect and the third-body perturbation up to the third order, i.e., including the harmonics in Equations (1) and (2) with \(l=2,3\). In order to investigate the long-term evolution, a doubly averaged model is used, i.e., the disturbing potential is averaged over the mean motion of the satellite and over the mean motion of the third body. The disturbing potential \(\mathcal{R}\) can be written as the sum of the contribution due to the second zonal harmonics, say \(\bar{\mathcal{R}}_{l_{2}}\), and the one due to Moon and Sun, say \(\bar{\mathcal{R}}_{\mathcal{Q}}\) and \(\bar{\mathcal{R}}_{\odot}\), respectively. We have  
\[\mathcal{R}=\bar{\mathcal{R}}_{l_{2}}+\bar{\mathcal{R}}_{\mathcal{Q}}\ +\bar{ \mathcal{R}}_{\odot}, \tag{3}\]  
where the first term in Equation (3) is the singly averaged secular oblateness effect [15]  
\[\bar{\mathcal{R}}_{l_{2}}=\frac{1}{4}I_{2}\frac{\mu_{\oplus}R_{\oplus}^{2}}{a ^{3}(1-e^{2})^{\frac{3}{2}}}\bigg{(}1-3\cos^{2}i\bigg{)}, \tag{4}\]where \(J_{2}\) is the second-order zonal coefficient of the Earth's gravity field harmonic expansion, \(\mu_{\odot}\) is the Earth's gravitational parameter and \(R_{\oplus}\) represents the equatorial mean radius of the Earth.