## 2 Theoretical Background
### Hamiltonian Dynamics  
In order to highlight the Hamiltonian structure of the problem, a coordinate change is required, by considering the Delaunay canonical variables \((L,G,H,\ell,g,h)\) defined as [16]  
\[ \begin{array}{l}L=\sqrt{\mu_{\oplus}a},\hskip 28.452756 ptG=L\sqrt{1-e^{2}}, \hskip 28.452756 ptH=G\cos i,\\ \ell=M,\hskip 56.905512 ptg=\omega,\hskip 56.905512 pth=\Omega.\end{array}  \tag{6}\]  
In this way, the Hamiltonian describing the long-term lunisolar effect on a Molinya satellite can be written as  
\[\mathcal{H}(L,G,H,\ell,g,h)=H_{kep}(L)+\mathcal{H}_{I_{2}}(G,H;L)+\mathcal{H}_ {\mathbb{Q}}\ (G,H,g,h;L)+\mathcal{H}_{\oplus}(G,H,g,h;L), \tag{7}\]  
where  
\[H_{kep}=-\frac{\mu_{\oplus}^{2}}{2L^{2}}, \tag{8}\]  
and  
\[\mathcal{H}_{I_{2}}=-\bar{\mathcal{R}}_{I_{2}},\hskip 14.226378 pt\mathcal{H}_{ \mathbb{Q}}\ =-\bar{\mathcal{R}}_{\mathbb{Q}}\,\hskip 14.226378 pt\mathcal{H}_{ \odot}=-\bar{\mathcal{R}}_{\odot}, \tag{9}\]  
written in terms of Delaunay variables. It has to be pointed out that in Equation (2), the harmonic argument vanishes for \(l=2\), \(p=1\) and \(m,s=0\). Since the corresponding term in \(\mathcal{H}_{\mathbb{Q}}\) only depends on the actions \((L,G,H)\), we will call this special harmonic the _lunar mean term_. As for the lunar case, for \(l=2\), \(p=1\) and \(m=0\), the solar harmonic argument in Equation (1) disappears and thus we will refer to the corresponding harmonic term as the _solar mean term_. Let us introduce the following notation:  
\[ \begin{array}{l}\mathcal{H}_{\mathbb{Q}}\ (G,H,g,h;L)=C_{0}^{\mathbb{Q}} \ \mathcal{A}_{0}^{\mathbb{Q}}\ (G,H;L)+\sum_{\alpha}C_{\alpha}^{\mathbb{Q}}\ \mathcal{A}_{\alpha}^{\mathbb{Q}}\ (G,H;L)\cos(\varphi_{ \alpha}^{\mathbb{Q}}\ ),\\ \mathcal{H}_{\odot}(G,H,g,h;L)=C_{0}^{\odot}\ \mathcal{A}_{0}^{\odot}(G,H;L)+ \sum_{\gamma}C_{\gamma}^{\odot}\mathcal{A}_{\gamma}^{\odot}(G,H;L)\cos( \varphi_{\gamma}^{\odot}),\end{array}  \tag{10}\]  
where:  
* \(\alpha\) and \(\gamma\) index the finite number of lunar and solar, respectively, harmonics retained in the model;
* \(C_{\alpha}^{\mathbb{Q}}\ \mathcal{A}_{\alpha}^{\mathbb{Q}}\ (G,H;L)\) is the \(\alpha\)-th lunar harmonic coefficient, and the constant term \(C_{\alpha}^{\mathbb{Q}}\) includes the lunar orbital parameters;
* \(C_{\gamma}^{\odot}\mathcal{A}_{\gamma}^{\odot}(G,H;L)\) is the \(\gamma\)-th solar harmonic coefficient and \(C_{\gamma}^{\odot}\) includes the solar orbital parameters;
* \(\alpha=0\) and \(\gamma=0\) denote the mean terms, i.e., \(C_{0}^{\mathbb{Q}}\ \mathcal{A}_{0}^{\mathbb{Q}}\ (G,H;L)\) is the _lunar mean term_ and \(C_{0}^{\odot}\mathcal{A}_{0}^{\odot}(G,H;L)\) is the _solar mean term_;
* \(\varphi_{\alpha}^{\mathbb{Q}}\) is the \(\alpha\)-th lunar argument and \(\varphi_{\gamma}^{\odot}\) is the \(\gamma\)-th solar argument. Generally, both sine and cosine trigonometric functions appear in the development of the lunar disturbing function. However, in our case, only cosine harmonics persist because of the values that \(s\) and \(y_{s}\) in Equation (2) assume when \(l=2,3\). The mean anomaly is cyclic in the doubly averaged Hamiltonian Equation (10); hence, the action \(L\) is a first integral, i.e., the semi-major axis is constant in the long term. Given the above considerations, the dynamics of the \(G,H,g\) and \(h\) are given by the following Hamilton equations:  
\[\left\{