%Version 3.1 December 2024
% See section 11 of the User Manual for version history
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%\documentclass[referee,sn-basic]{sn-jnl}% referee option is meant for double line spacing

%%=======================================================%%
%% to print line numbers in the margin use lineno option %%
%%=======================================================%%

%%\documentclass[lineno,pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%=========================================================================================%%
%% the documentclass is set to pdflatex as default. You can delete it if not appropriate.  %%
%%=========================================================================================%%

%%\documentclass[sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%Note: the following reference styles support Namedate and Numbered referencing. By default the style follows the most common style. To switch between the options you can add or remove Numbered in the optional parenthesis. 
%%The option is available for: sn-basic.bst, sn-chicago.bst%  
 
%%\documentclass[pdflatex,sn-nature]{sn-jnl}% Style for submissions to Nature Portfolio journals
%%\documentclass[pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style
%\documentclass[pdflatex,sn-mathphys-num]{sn-jnl}% Math and Physical Sciences Numbered Reference Style
%%\documentclass[pdflatex,sn-mathphys-ay]{sn-jnl}% Math and Physical Sciences Author Year Reference Style
%%\documentclass[pdflatex,sn-aps]{sn-jnl}% American Physical Society (APS) Reference Style
%%\documentclass[pdflatex,sn-vancouver-num]{sn-jnl}% Vancouver Numbered Reference Style
\documentclass[pdflatex,sn-vancouver-ay]{sn-jnl}% Vancouver Author Year Reference Style
%%\documentclass[pdflatex,sn-apa]{sn-jnl}% APA Reference Style
%%\documentclass[pdflatex,sn-chicago]{sn-jnl}% Chicago-based Humanities Reference Style


%\documentclass[sn-authoryear]{sn-jnl} %Este es viejo.


%%%% Standard Packages
%%<additional latex packages if required can be included here>

\usepackage{graphicx}%
\usepackage{multirow}%
\usepackage{amsmath,amssymb,amsfonts}%
\usepackage{amsthm}%
\usepackage{mathrsfs}%
\usepackage[title]{appendix}%
\usepackage{xcolor}%
\usepackage{textcomp}%
\usepackage{manyfoot}%
\usepackage{booktabs}%
\usepackage{algorithm}%
\usepackage{algorithmicx}%
\usepackage{algpseudocode}%
\usepackage{listings}%
%%%%


\usepackage{ulem} % Necesario para \sout
\usepackage{natbib} % Para la bibliografia con apellidos y años
\usepackage{wasysym} %Primero amsmath. Si no problemas.
\usepackage{xcolor}
%\usepackage[toc,page]{appendix}
%%%%
\newcommand{\sign}{\mathop{\rm sign}\nolimits}
\newcommand{\Eqref}[1]{(\ref{#1})}
\newcommand{\Frac}[2]{\mbox{$\displaystyle\frac{#1}{#2}$}}
\newcommand{\Int}{\displaystyle\int}
\newcommand{\Sqrt}[1]{\displaystyle\sqrt{#1}}


%%%%%=============================================================================%%%%
%%%%  Remarks: This template is provided to aid authors with the preparation
%%%%  of original research articles intended for submission to journals published 
%%%%  by Springer Nature. The guidance has been prepared in partnership with 
%%%%  production teams to conform to Springer Nature technical requirements. 
%%%%  Editorial and presentation requirements differ among journal portfolios and 
%%%%  research disciplines. You may find sections in this template are irrelevant 
%%%%  to your work and are empowered to omit any such section if allowed by the 
%%%%  journal you intend to submit to. The submission guidelines and policies 
%%%%  of the journal take precedence. A detailed User Manual is available in the 
%%%%  template package for technical guidance.
%%%%%=============================================================================%%%%

%% as per the requirement new theorem styles can be included as shown below
\theoremstyle{thmstyleone}%
\newtheorem{theorem}{Theorem}%  meant for continuous numbers
%%\newtheorem{theorem}{Theorem}[section]% meant for sectionwise numbers
%% optional argument [theorem] produces theorem numbering sequence instead of independent numbers for Proposition
\newtheorem{proposition}[theorem]{Proposition}% 
%%\newtheorem{proposition}{Proposition}% to get separate numbers for theorem and proposition etc.

\theoremstyle{thmstyletwo}%
\newtheorem{example}{Example}%
\newtheorem{remark}{Remark}%

\theoremstyle{thmstylethree}%
\newtheorem{definition}{Definition}%

\raggedbottom
%%\unnumbered% uncomment this for unnumbered level heads

\begin{document}
	
\title[\textcolor{red}{Quasi-frozen polar orbits and resonance effects in low-altitude orbits around the Moon}]{Quasi-frozen polar orbits and resonance effects in low-altitude orbits around the Moon}


%%=============================================================%%
%% GivenName	-> \fnm{Joergen W.}
%% Particle	-> \spfx{van der} -> surname prefix
%% FamilyName	-> \sur{Ploeg}
%% Suffix	-> \sfx{IV}
%% \author*[1,2]{\fnm{Joergen W.} \spfx{van der} \sur{Ploeg} 
%%  \sfx{IV}}\email{iauthor@gmail.com}
%%=============================================================%%


\author*[1]{\sur{Jean Paulo S. Carvalho}}\email{jean.carvalho@ufrb.edu.br}

\author[2,3]{\sur{Daniel Casanova}}\email{casanov@unizar.es}
\equalcont{These authors contributed equally to this work.}

\affil*[1] {\orgdiv{Centro de Ci\^encia e Tecnologia em Energia e Sustentabilidade}, \orgname{Universidade Federal do Rec\^oncavo da Bahia}, \orgaddress{\street{Av. Centen\'ario}, \city{Feira de Santana}, \postcode{44042280}, \state{Bahia}, \country{Brazil}}} 

\affil[2]{\orgdiv{Centro Universitario de la Defensa de Zaragoza}, \orgname{Academia General Militar}, \orgaddress{\street{Crta. de Huesca s/n}, \city{Zaragoza}, \postcode{50090}, \state{Arag\'on}, \country{Spain}}}

\affil[3]{\orgdiv{IUMA}, \orgname{Universidad de Zaragoza}, \orgaddress{\street{Pedro Cerbuna, 12}, \city{Zaragoza}, \postcode{50009}, \state{Arag\'on}, \country{Spain}}}


%%==================================%%
%% Sample for unstructured abstract %%
%%==================================%%

\abstract{In this work, we analyze orbits around the Moon, considering the primary perturbing forces, including the effects of the selenopotential, third-body perturbation and solar radiation pressure (SRP). We adopt a simplified model with zonal truncation up to order eighteen, incorporating sectoral and tesseral terms of order and degree three, named Model 18 $\times$ 3. Additionally, we develop an analytical model that includes zonal harmonics up to the $J_{18}$ term, named Model 18 $\times$ 0. With the analytical model, we analyze the conditions of frozen orbits when the perilune argument is 270 and 90 degrees. We show the results in graphs of eccentricity versus inclination, where we obtain initial conditions of frozen orbits. Numerically integrating the equations of motion using the initial conditions obtained in the analytical model, we identify quasi-frozen polar orbits and regions of large eccentricity variation that can aid in the removal of space debris. To validate the results, we used the GMAT software considering the selenopotential model of order and degree 50 $\times$ 50.}

%%================================%%
%% Sample for structured abstract %%
%%================================%%

% \abstract{\textbf{Purpose:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Methods:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Results:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Conclusion:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.}

\keywords{Lunar orbit; Selenopotential; Frozen orbit; Resonance}

%%\pacs[JEL Classification]{D8, H51}

%%\pacs[MSC Classification]{35A01, 65L10, 65L12, 65L20, 65L70}

\maketitle

\section{Introduction}\label{sec1}

Currently, several nations are planning to return to the Moon to establish a lunar base and/or further explore this celestial body, as exemplified by the Artemis Program. Notably, one of the mission's key objectives is to land the first woman on the Moon in 2026. Several space vehicles have already orbited the Moon, and new missions are underway. For example, Firefly's Blue Ghost Mission 1, named Ghost Riders in the Sky, launched on January 15, 2025, and successfully completed its 45-day journey to the Moon spending 25 days orbiting Earth, 4 days in transit, and 16 days in lunar orbit before making a soft landing on March 2, 2025. Blue Ghost will then conduct 14 days of surface operations. In this study, we analyze regions around the Moon that could support long- and short-lifetime orbits. Regions with short-lifetime orbits may serve as potential solutions for space debris removal.

For this analysis, we consider the main perturbing forces in the dynamical model, namely, the non-uniform mass distribution of the Moon, solar radiation pressure (SRP), and third-body perturbations from the Earth. In \citet{Oesterwinter}, the author develops a semi-analytical theory of lunar satellite motion using the Von Zeipel perturbation method, analyzing secular and long-periodic terms up to the second order. The primary perturbations considered include the Moon's non-spherical gravitational field and the gravitational attraction of the Earth and the Sun. Similarly, in \citet{Giacaglia}, the authors apply the Von Zeipel method to develop a semi-analytical theory of lunar satellite motion, incorporating perturbations from the Moon's non-spherical gravity, Earth-Moon-Sun gravitational interactions, SRP, and the Moon's physical libration. Additionally, in \citet{Carvalho6}, using the Lie-Hori method, the disturbing potential due to the Moon's gravitational non-sphericity is derived up to the second order. Moreover, in \citet{Saedeleer,Henrard}, an analytical theory for artificial lunar satellites is presented, employing the Lie Transform perturbation method to compute the averaged Hamiltonian of the problem. In \citet{Ren}, the authors derive an analytical solution for the three-body perturbation dynamics using the double-averaged model, considering a perturbing body in an inclined elliptical orbit, and in \citet{Nie2026} with high elliptical orbits. In \citet{Abad}, analytical theories based on the Lie-Deprit transform are used to determine families of periodic orbits for a lunar orbiter, applying the single averaging method to derive the equations of motion. In \citet{Sharma}, different lunar potential models are compared, incorporating high-order spherical harmonics (zonal and tesseral terms). The study examines the effect of the right ascension of the ascending node on orbital lifetime variations with altitude, particularly for polar and circular lunar orbits, using the lunar orbit propagation software ``lprop.'' Additionally, the lifetime of high-altitude circular lunar orbits is analyzed. In \citet{Liu}, a low-thrust process is employed to implement an orbit maintenance strategy for ultra-low-altitude, low-inclination lunar orbits. In \citet{Tzirti}, the secular dynamics of lunar satellites are investigated using high-precision gravitational models and frequency analysis techniques based on the Prony method.  
In \citet{Pontani}, the orbital dynamics of a primary microsatellite and its associated nanosatellites post-launch are analyzed using non-singular equinoctial elements for orbit propagation and a numerical averaging technique. The dynamical model incorporates the lunar gravitational field LPE100K, considering all harmonics up to order 70, as well as third-body perturbations from the Earth and the Sun. In  \citet{Folta}, the authors performed a frozen orbit analysis considering two orbital regimes: low lunar orbits and elliptical orbits. The model includes third-body gravitational perturbations from the Earth, the Sun, and Jupiter, as well as accelerations due to the lunar gravitational potential represented through spherical harmonics up to degree and order 165 or 100, using the LP165 and LP100K4 models, respectively. In \citet{Konopliv}, a lunar gravity field model up to degree and order 60 is presented, developed using two-way and three-way Doppler tracking data obtained from the Lunar Orbiter I to V missions.  

Numerical studies on the characteristics of lunar orbits are presented in works such as \citet{Popandopulo, Sharma, Melo, Liana}. Research on frozen orbits for artificial lunar satellites includes studies such as \citet{Abad, Elipe, Carvalho1, Khattab, Magdy, Nie, ElSalam, Xu2025}. In \citet{Gomes}, the lifetime of orbits around natural satellites, including the Moon, that move in elliptical trajectories around their host planets, is analyzed. The results are presented as maps illustrating orbital lifetime as a function of initial inclination and eccentricity. In \citet{Magdy}, moderate-altitude frozen orbits around the Moon are investigated using an averaged Hamiltonian, considering zonal harmonics up to order 7 (7x0) and third-body perturbations from the Earth, modeled as an inclined elliptical orbit.  

This work is organized as follows: in Section \ref{sec2}, we introduce the mathematical model of the lunar potential, which incorporates zonal coefficients up to the eighteen order, as well as sectoral and tesseral terms up to the third order and degree. We then present the disturbing potential due to gravitational perturbations from the Earth, based on the single averaged model in an inclined elliptical orbit, as presented in \citet{Carvalho3}. Additionally, we adopt the formulation for the SRP model as presented in \citet{Carvalho} and \citet{Carvalho7}. In Section \ref{res}, we examine the conditions of frozen polar orbits and the resonance effect due to SRP coupling and the oblateness of the Moon. Here, we consider the analytical model to find the initial conditions of frozen orbits and the numerical model to propagate these orbits and analyze their long-term evolution. We perform comparisons of the numerical model with the GMAT software. In the final section, we summarize our main findings and discuss their implications for lunar mission planning and the sustainability of the Moon's orbital region.




\section{Dynamical model}\label{sec2}

In this work, a lunar orbiter is being considered. Consequently, the main perturbative effects are caused by the Moon, acting as the central body, the Earth as third body, and the Solar Radiation Pressure (SRP). Thus, the problem's potential is described as follows:

\begin{equation}\label{0}
	\begin{array}{l}
		\displaystyle
		U = U_{\leftmoon} + U_{\oplus} + U_{SRP},
	\end{array}
\end{equation}
where $U_{\leftmoon}$ denotes the gravitational potential of the Moon, $U_{\oplus}$ correspond to the gravitational potential exerted by the Earth. Here, $U_{SRP}$ represents the potential caused by the Solar Radiation Pressure. It is worth noting that $U_{SRP}$ is a non-conservative force; however, under certain conditions, it can be treated as an effectively conservative force, allowing its potential to be expressed.

We present the influence of different perturbations on the satellite/debris in  Figure~\ref{fig:AccVSDis}. It shows the acceleration caused by the described perturbations plus others, depending on the altitude of the satellite/debris with respect to the Moon's surface, and justify the selection of our simplified model.  

\begin{figure}[!h]
	\centering
	\resizebox{0.99\columnwidth}{!}{\includegraphics{accVSdistance.pdf} }
	\caption{Acceleration (km/s$^2$) caused by the different perturbations depending on the altitude of the satellite/debris with respect to the Moon surface (km).}
	\label{fig:AccVSDis}       % Give a unique label
\end{figure}

\subsection{Disturbing potential of the central body (the Moon) and its averaged formulation}

The problem of a satellite/debris orbiting the Moon is schematically presented in Figure~\ref{fig:GeometryMoon}, where the gross geometry of the problem is shown. Since the central body rotates, it becomes necessary to define a direction which defines this rotation with respect to the inertial system. Therefore the line $ O-Md $ defines the instantaneous direction of some specific meridian, therefore it rotates according to the spin rate of the central body, which is named $\nu$.  The longitude $\lambda_M$ is measured in the equatorial plane from the rotating meridian ($ O-Md $) to the spacecraft meridian, in the equator plane. The ascending node is defined relative to the inertial $OX$ axis, and its longitude is denoted by $h$. The classical orbital elements are: $a$ (semi-major axis), $e$ (eccentricity), and $i$ (inclination). The argument of periselene is $g$, the true anomaly is $f$, and the argument $f + g$ defines the angular position of the spacecraft in the orbital plane. The angle $\phi$ represents the latitude of the S/C.

%
\begin{figure}[!htb]
	\centering
	\resizebox{0.99\columnwidth}{!}{\includegraphics{GeometryMoon.pdf} } 
	\caption{Reference geometry for a S/C orbiting the Moon. Adapted from~\citet{Carvalho3}}
	\label{fig:GeometryMoon}       % Give a unique label
\end{figure}
%
From Figure~\ref{fig:GeometryMoon}, the following relations can be derived, which are used in subsequent computations:
\begin{eqnarray}
	\nonumber \lambda_M  =  x+y, \\
	\nonumber h = x + \nu t , \\
	\cos(f+g)  =  \cos y \cdot \cos\phi, \\
	\nonumber \sin (\phi) =\sin i \cdot \sin(f+g), \\
	\nonumber \sin (f+g) \cdot \cos i = \sin y \cdot \cos\phi.
\end{eqnarray}

Due to rotational motion, in general, celestial bodies are not a perfect sphere, but have the shape of an oblate spheroid with an equatorial diameter that exceeds the polar diameter. The resulting equatorial bulge exerts a force that pulls the artificial satellite back to the equatorial plane whenever it is above or below this plane and thus tries to align the orbital plane with the equator \citep{Montenbruck}. In this work, the Moon is considered as the celestial body acting as the central object and, to mathematically describe this disturbing effect, expansions in spherical harmonics via Legendre polynomials are used. As described in \citet{Montenbruck}, the gravitational potential is given by:

\begin{equation}\label{1}
	\begin{array}{l}
		\displaystyle
		U_{\leftmoon}  =\frac{\mu}{r}\sum_{n=2}^{\infty}\sum_{m=0}^{n} \left( \frac{R_{M}}{r} \right)^{n} P_{nm}(\sin\phi) \, (C_{nm} \, \cos(m\lambda_M)+S_{nm} \, \sin(m\lambda_M)).
	\end{array}
\end{equation}

Here, $\mu = 4902.801056 \,$ km$^3 / $s$^2$ is the gravitational constant of the Moon, $R_{M} = 1738$ km is the equatorial radius of the Moon, $r=||\mathbf{r}||$, where $\mathbf{r}$ is the position vector of the satellite/debris, the angle $\phi$ is the latitude of the orbit with respect to the equator of the Moon and the angle $\lambda_M$ is the longitude (counted positively towards the East) of the point $r$ with respect to the line  $O-Md$. Finally, $P_{nm}$ is the  associated Legendre polynomial of degree $n$ and order $m$, $C_{nm}$ and $S_{nm}$ are the coefficients that describe the dependence of the internal distribution of the Moon mass (spherical harmonic coefficients). In particular, $J_{n}=-C_{n0}$. In the case of the Moon, the selenopotential presents a reduction in the values of the coefficients $C_{nm}$ and $S_{nm}$  as the order and degree of the polynomials and associated functions increases, causing rounding errors when calculated on the computer, due to the limited number of significant digits. To minimize these effects, it is customary to use the coefficients in a completely normalized form \citet{Kuga}. The normalized coefficients $\overline{C}_{nm}$ and $\overline{S}_{nm}$ are defined as \citet{Kuga,Vinti,Roncoli}:


\begin{equation}
	%\begin{array}{l}
	\binom{\overline{C}_{nm}}{\overline{S}_{nm}}=\sqrt{\frac{(n+m)!}{(4n+2)(n-m)!}}\binom{{C}_{nm}}{{S}_{nm}}.
	%\end{array}
	\label{2}
\end{equation}

Making use of the normalized coefficients, the Selenopotential can be rewritten as~\citep{Kuga,Vinti}:

\begin{equation}\label{4}
	\begin{array}{l}
		\displaystyle
		U_{\leftmoon}  =\frac{\mu}{r}\sum_{n=0}^{\infty}\sum_{m=0}^{n}
		\left( \frac{R_M}{r} \right)^n \overline{P}_{nm}(\sin\phi)(\overline{C}_{nm} \cos(m\lambda_{M})+
		\overline{S}_{nm} \sin(m\lambda_{M})),
	\end{array}
\end{equation}
where

\begin{equation}\label{5}
	%\begin{array}{l}
	\overline{P}_{nm}=\sqrt{\frac{(4n+2)(n-m)!}{(n+m)!}} \, \, {P}_{nm}.
	%\end{array}
\end{equation}


In this work, we analyze the effect of the main terms of spherical harmonics according to the values presented in \citet{Roncoli}, see pages 59 and 60.
Table ~\ref{tab:LP150Q} presents the values of the spherical harmonics considered, however the zonal terms are presented up to the $J_{9}$ term, the other zonal harmonics can be consulted in \citet{Roncoli}.

% For tables use

\begin{table}
	\centering
	% table caption is above the table
	\caption{LP150Q gravitational field, \cite{Roncoli}. Where $\mu=4902.801076$ km$^{3}/$s$^{2}$ and $R_{M}=1738$ km,
		$\rho=2.649072347 10^{-6}$ rad/s (Moon rotation rate).}
	\label{tab:LP150Q}       % Give a unique label
	% For LaTeX tables use
	\begin{tabular}{c|rr}
		\hline\noalign{\smallskip}
		Coefficient & Normalized  \\
		\noalign{\smallskip}\hline\noalign{\smallskip}
		$J_{2}$ & $0.909010949481 \times 10^{-4}$ \\
		$J_{3}$ & $0.320307167959 \times 10^{-5}$  \\
		$J_{4}$ & $-0.321409545028 \times 10^{-5}$  \\
		$J_{5}$ & $0.221009876393 \times 10^{-6}$  \\
		$J_{6}$ & $-0.376479064475 \times 10^{-5}$  \\
		$J_{7}$ & $-0.561330656403 \times 10^{-5}$  \\
		$J_{8}$ & $-0.231953905520 \times 10^{-5}$  \\
		$J_{9}$ & $0.354241582136 \times 10^{-5}$  \\
		$C_{22}$ & $0.346376274208 \times 10^{-4}$  \\
		$C_{31}$ & $0.263418358622 \times 10^{-4}$  \\
		$C_{32}$ & $0.141853316786 \times 10^{-4}$  \\
		$C_{33}$ & $0.122862645044 \times 10^{-4}$  \\
		$S_{22}$ & $0.144063503540 \times 10^{-7}$  \\
		$S_{31}$ & $0.546307860882 \times 10^{-5}$  \\
		$S_{32}$ & $0.488913911795 \times 10^{-5}$  \\
		$S_{33}$ & $-0.178246270720 \times 10^{-5}$  \\
		\noalign{\smallskip}\hline
	\end{tabular}
\end{table}

Now, the disturbing potential due to the non-uniform distribution of mass of the central body given by Eq.~(\ref{4}) is rewritten as a function of the orbital elements considering the zonal harmonics $J_{2}-J_{18}$ and the terms $C_{22}$, $C_{31}, C_{32}, C_{33}$, $S_{22}, S_{31}, S_{32}, S_{33}$ (Note that $C_{21}=S_{21} =0$, and $S_{20}=S_{30}=0$). Later, a single averaged method is applied to eliminate the short-periodic terms. See \citet{Carvalho3} to show the development of the equations. However, this work explicitly presents the expansion of the zonal coefficients up to the $J_{9}$, while in \citet{Carvalho3} it was presented up to the $J_{6}$. It should be noted that the normalization factors calculated by Eq.~(\ref{5}) have now also been included. Note that, the equatorial plane of the Moon is considered as the reference plane, and consequently, the averaged disturbing potential due to the non-uniform distribution of mass of the central body is given by:


\[
\langle U_{\leftmoon} \rangle = \langle R_{J2} \rangle +  \langle R_{J3} \rangle + \cdots + \langle R_{J18} \rangle +  \langle R_{C22,S22} \rangle  + 	\langle R_{C31,S31} \rangle +  \langle R_{C32,S32} \rangle + \langle R_{C33,S33} \rangle.
\]
See Appendix \ref{secA1} for a development of terms up to ninth order (Model 9 $\times$ 3).


\subsection{Perturbation due to third body and its averaged formulation}

In the perturbation caused by the third body it is employed the classical expansion of Legendre polynomials, retaining only second order terms in the distance ratio. However, in the derivation of the disturbing potential, it is used the single averaged model, wherein the third body is situated in an elliptical and inclined orbit. The gravitational influence of the Earth is expressed through the perturbation of the third body within the inertial frame Oxyz centered on the Moon. Thus, it is obtained:


\begin{equation}\label{3}
	\begin{array}{l}
		U_{3b}= -\mu_{d}\left(\frac{\mathbf{r}-\mathbf{r}_{d}}{\|\mathbf{r}-\mathbf{r}_{d}\|^3}+
		\frac{\mathbf{r}_{d}}{\|\mathbf{r}_{d}\|^3}\right),
	\end{array}
\end{equation}

where $\mu_{d} = G m_{d}$ represents the gravitational parameter of the disturbing body. The vector $\mathbf{r_{d}}$ is the position vector of the disturbing body with respect to central body and $\mathbf{r}$ represents the position of the satellite/debris. The averaged perturbation due to third body is $ \langle U_{3b} \rangle $, where the term $3b$ could be the Sun or the Earth (See \citet{Carvalho, Carvalho3} for a detailed explanaition).


\subsection{Perturbation due to Solar Radiation Pressure and its averaged formulation}

In this work, it is considered the solar radiation pressure (SRP) on a lunar orbit, as other authors did in previous works \citet{Kubo-oka, Kubo}. The model presented by \citet{Kubo-oka} assumed   the satellite as a regular octagonal cylinder. However, in \citet{Kubo} the authors presented an enhanced model of the SRP of the relay satellite, where the shape of the relay satellite is more realistic. More precisely, it is modeled as a combination of a regular octagonal cylinder corresponding to the main body of the spacecraft and a column corresponding to the mounting assembly with the lunar orbiter. In the SRP model we assume that the artificial satellite remains in a fixed orientation perpendicular to the Sun at all times. This satellite can be coupled to a solar sail to increase the area-to-mass ratio, thereby enhancing the effect of the SRP \citep{Tresaco,Tresacob,McInnes}. One possible application could be the removal of space debris from lunar orbit. Note that the force on a solar sail depends on the sail's area and orientation. When the perturbation term due to solar radiation pressure is added to the system, it only maintains its Hamiltonian character when the sail is aligned with respect to the Sun (i.e. no sail effect) or in the case where the sail is perpendicular to the Sun-sail direction (i.e. maximum sail effect). For other navigation directions, the systems is no longer Hamiltonian. For reasons of convenience, we neglect the shading effect, as it was also neglected in \citet{Kubo-oka, Kubo}. In \citet{Carvalho}, a single averaged model of the SRP is presented, considering the Sun in an elliptical and inclined orbit, based on the approach from \citet{Tresacob,Tresaco}. To avoid redundancy, we present only a portion of the development here, with further details available in \citet{Carvalho}. As discussed in \citet{Tresacob,Tresaco}, given these considerations, the SRP can be expressed as a perturbing potential, neglecting the Earth's shadow, as follows:
\begin{equation}\label{U_SRP}
	\begin{array}{l}
		U_{SRP}=-\beta\mu_{\odot}\left(\frac{1}{\|\mathbf{r}-\mathbf{r_{\odot}}\|}\right),
	\end{array}
\end{equation}
where $\beta$ (see \citet{McInnes,Tresacob,Tresaco} for more details) is defined as:

\begin{equation}\label{betaEq}
	\begin{array}{l}
		\beta= 1.53 \cdot \frac{A}{m},
	\end{array}
\end{equation}

where the value 1.53 $ \text{gr/m}^2 $ is a factor computed considering the Sun's luminosity and the speed of light for the case of a solar sail around the Earth/Moon~\citep{McInnes}. Here, $A$ is the sail area and $m$ is the total mass of the spacecraft, $\mu_{\odot}$ represents the gravitational parameter of the Sun, and $\mathbf{r_{\odot}}$ is the position vector of the Sun with respect to Earth-Moon.

Thus, taking into account that the sail is assumed to maintain an orientation perpendicular to the Sun line and is located at a distance from the Sun $||\mathbf{r}-\mathbf{r}{\odot}||$, and that the case of interest is $r \ll r{\odot}$, the following term is expanded as follows:

\begin{equation}\label{10}
	\begin{array}{ll}
		\displaystyle
		\frac{1}{\|\mathbf{r}-\mathbf{r}_{\odot}\|} & = \displaystyle \frac{1}{r_{\odot}}\left(1-2\frac{r}{r_{\odot}}\cos\psi+
		\left(\frac{r}{r_{\odot}}\right)^2\right)^{-1/2}\\[1.5ex]
		& = \displaystyle \frac{1}{r_{\odot}}\sum_{n\geq0}^2\left(\frac{r}{r_{\odot}}\right)^nP_n(\cos\psi).
	\end{array}
\end{equation}
Then, the potential function can be expanded in terms of Legendre polynomials up to the second order:

\begin{equation}\label{11}
	\begin{array}{l}
		\displaystyle
		U_{SRP}=-\beta \, \frac{\mu_{\odot}}{r_{\odot}}\left[1+\frac{r}{r_{\odot}}\cos\psi+
		\frac{1}{2} \left(\frac{r}{r_{\odot}}\right)^2 (3\cos^2\psi-1)  \right],
	\end{array}
\end{equation}
where $\psi$ is the angle between the vectors $\mathbf{r}$ and $\mathbf{r_{\odot}}$, which can be expressed in terms of the orbital elements of the Sun and the spacecraft. As previously commented, a more detailed explanation is presented in \citet{Carvalho}. It is assumed that the apparent motion of the Sun follows an elliptical and inclined orbit around the Earth-Moon with semimajor axis $a_{\odot}$, eccentricity $e_{\odot}$, inclination $i_{\odot}$, periselene argument $g_{\odot}$, longitude of the ascending node $h_{\odot}$, the mean motion $n_{\odot}$ and $\lambda_{\odot}$ the longitude of the Sun in the ecliptic plane. However, to carry out an analytical study, it is considered $e_{\odot}=0$ %in equation (22) of \citet{Carvalho} 
and it is obtained the following averaged SRP potential:


\begin{equation}\label{12}
	\begin{array}{rl}
		\langle U_{SRP}\rangle & ={\displaystyle\frac {3}{8}}e {n^2_{{\odot}}}\beta\,a a_{{\odot}}
		\Big(
		(1+\cos i)( 1+\cos \epsilon) \cos(g+h-\lambda_{\odot}) \\[1.5ex] %Seria \psi_1
		&-( 1+\cos \epsilon ) ( \cos i -1) \cos( g-h+\lambda_{\odot}) \\[1.5ex] %Sería -\psi_2
		&+ 2\, \cos( g-\lambda_{\odot}) \sin  i  \sin  \epsilon   - %Sería \psi_3
		2\, \cos( g+\lambda_{\odot}) \sin  i  \sin  \epsilon  \\[1.5ex] %Sería \psi_4
		&-(\cos\epsilon -1)( 1+\cos i) \cos(g+h+\lambda_{\odot}) \\[1.5ex] %Sería \psi_5
		&+ (\cos\epsilon -1 )( \cos i -1) \cos( g-h-\lambda_{\odot}) \Big), %Sería -\psi_6
		\Big.
	\end{array}
\end{equation}
where the following resonances appears:

\begin{equation}\label{13}
	\begin{array}{llll}
		\psi_1 = h+g-\lambda_{\odot}, & \qquad \psi_2 = h-g-\lambda_{\odot}, & \qquad \psi_3 = g-\lambda_{\odot},&\\
		\psi_4 = g+\lambda_{\odot}, & \qquad \psi_5 = h+g+\lambda_{\odot}, & \qquad \psi_6 = h-g+\lambda_{\odot}.&
	\end{array}
\end{equation}
Thus, the averaged disturbing potential due to the SRP is expressed by considering the six resonant angles mentioned previously:
\begin{equation}\label{13b}
	\begin{array}{rl}
		\langle U_{SRP}\rangle & ={\displaystyle\frac {3}{8}}e {n^2_{{\odot}}}\beta\,a a_{{\odot}}
		\Big( (1+\cos i)( 1+\cos \epsilon) \cos(\psi_1) %Seria \psi_1
		-( 1+\cos \epsilon ) ( \cos i -1) \cos(-\psi_2)  \\ %[1.5ex] %Sería -\psi_2
		& + 2\, \cos(\psi_3) \sin  i  \sin  \epsilon   - %Sería \psi_3
		2\, \cos(\psi_4) \sin  i  \sin  \epsilon -  %Sería \psi_4
		(\cos\epsilon -1)( 1+\cos i) \cos(\psi_5) \\ %[1.5ex] %Sería \psi_5
		& + (\cos\epsilon -1 )( \cos i -1) \cos(-\psi_6) %Sería -\psi_6 
		\Big).
	\end{array}
\end{equation}

Consequently, it is possible to decompose the average disturbing potential caused by SRP into six distinct resonant terms:

\begin{equation}\label{eq:resonant_equations}
	\begin{array}{lrl}
		\langle R_{SRP_{1}} \rangle =&  \mathcal{K} &
		\left(1+\,\cos \left( i \right)  \right)  \left( 1+\cos \left( \epsilon \right)  \right) \cos \left( \psi_1 \right), \\
		\langle R_{SRP_{2}} \rangle =&  - \mathcal{K} &
		\left( 1\,\cos \left( i \right) -1 \right) \left( 1+\cos \left( \epsilon \right)  \right)
		\cos \left( \psi_2 \right), \\
		\langle R_{SRP_{3}} \rangle =&  2 \mathcal{K} &
		\cos \left( \psi_3 \right) \sin \left( i \right) \sin \left( \epsilon \right), \\
		\langle R_{SRP_{4}} \rangle =& - 2  \mathcal{K} &
		\cos \left(\psi_4 \right) \sin \left( i \right) \sin  \left( \epsilon \right), \\
		\langle R_{SRP_{5}} \rangle =&  -  \mathcal{K} &
		\left( \,\cos \left( \epsilon \right) -1 \right)  \left( 1+\cos \left( i \right)  \right) \cos \left(\psi_5 \right), \\
		\langle R_{SRP_{6}} \rangle =& \mathcal{K} &
		\left( \,\cos \left( i \right) -1 \right)  \left( \cos \left( \epsilon \right) -1 \right) \cos \left(\psi_6\right),
	\end{array}
\end{equation}
where $\mathcal{K} = \frac {3}{8}\,e \,{n_{\odot}}^{2} \beta \, a \, a_{\odot},$ and $\epsilon= i_{\odot}$.








\section{Frozen and resonant orbits}\label{res}
This section examines the conditions required to achieve frozen and resonant orbits, with emphasis on the key factors governing their stability and long-term behavior. Frozen orbits are characterized by orbital elements that remain nearly constant over time, whereas resonant orbits arise from specific commensurabilities between the orbital periods of two or more bodies. A proper understanding of these orbital configurations is essential for the analysis of long-term dynamical behavior.

The General Mission Analysis Tool (GMAT) is an open-source software package developed by NASA for the simulation and analysis of space missions. It is widely used for tasks such as orbital trajectory design, maneuver planning, flight optimization, and space systems modeling. GMAT supports mission analysis across a broad range of flight regimes, from low Earth orbit to lunar, libration point, and deep-space missions. Designed to support real missions, including the Lunar Reconnaissance Orbiter, GMAT combines high accuracy, flexibility, and extensibility. In this work, GMAT was used with a lunar gravity field (selenopotential) model of degree and order 50, providing a reference for comparison with our simplified dynamical model.


\subsection{Analytical model: frozen orbits}
By considering the zonal harmonic terms $J_{2} + \cdots + J_{18}$ in the perturbing potential---which are the dominant contributors to the spacecraft dynamics (see Table~\ref{tab:LP150Q}, where terms up to $J_{9}$ are listed; higher-order coefficients can be found in \cite{Roncoli}, p.~59)---and fixing the argument of periapsis $g$ at either $90$ degrees ($\frac{\pi}{2}$ radians) or $270$ degrees ($\frac{3\pi}{2}$ radians), the Lagrange Planetary Equations (LPE), specifically $\dot{e}$ and $\dot{i}$, vanish, while $\dot{g} \neq 0$.

As is well known, frozen orbits require the simultaneous satisfaction of the conditions $\dot{e}=0$, $\dot{i}=0$, and $\dot{g}=0$. Consequently, we impose $\dot{g}=0$ and analyze the resulting relationship between eccentricity and inclination, $e=f(i)$, using contour plots, as shown in Figures ~\ref{fig:Fig599}, ~\ref{fig:Fig601}, ~\ref{fig:Fig602}, ~\ref{fig:Fig603}, ~\ref{fig:Fig604} and ~\ref{fig:Fig605}. The blue and red contours indicate zero-level curves, which correspond to frozen-orbit solutions (see, for example, \cite{Abad} and \cite{Kikuchi}). 

To generate these figures, we included the zonal harmonics of the perturbing potential as well as the third-body perturbation. To simplify the dynamical model, the third body was represented using a double-averaged formulation assuming an elliptical and planar orbit. It is worth noting that \cite{Abad} adopted a circular and planar model, whereas \cite{Sirwah} and \cite{Nie2024} considered an elliptical and inclined configuration.

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig599.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=1838$~km. The model includes the zonal harmonics $J_{2}+ J_{7}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit. The blue and red curves correspond to $g=\pi/2$ and $g=3\pi/2$, respectively.}
		\label{fig:Fig599}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig601.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=2238$~km. The model includes the zonal harmonics $J_{2} + J_{7}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit. The blue and red curves correspond to $g=\pi/2$ and $g=3\pi/2$, respectively.}
		\label{fig:Fig601}
	\end{minipage}
\end{figure}

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig602.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=2238$~km. The model includes the zonal harmonics $J_{2}+ \cdots + J_{7}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit.}
		\label{fig:Fig602}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig603.eps}
		\caption{Zoomed view of Fig.~\ref{fig:Fig602}, highlighting the frozen-orbit solutions. The blue and red curves correspond to $g=\pi/2$ and $g=3\pi/2$, respectively.}
		\label{fig:Fig603}
	\end{minipage}
\end{figure}

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig604.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=2238$~km. The model includes the zonal harmonics $J_{2}+ \cdots + J_{18}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit.}
		\label{fig:Fig604}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig605.eps}
		\caption{Zoomed view of Fig.~\ref{fig:Fig604}, highlighting the frozen-orbit solutions. The blue and red curves correspond to $g=\pi/2$ and $g=3\pi/2$, respectively.}
		\label{fig:Fig605}
	\end{minipage}
\end{figure}

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig606.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=1838$~km. The model includes the zonal harmonics $J_{2}+ \cdots + J_{18}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit.}
		\label{fig:Fig606}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig608.eps}
		\caption{Frozen-orbit conditions for a lunar orbiter with semi-major axis $a=1938$~km. The model includes the zonal harmonics $J_{2}+ \cdots + J_{18}$ and the third-body disturbing potential (Earth), assuming an elliptical and planar orbit.}
		\label{fig:Fig608}
	\end{minipage}
\end{figure}

To validate our model, we compared our results with those reported in previous works. Considering the $J_{2}+J_{7}$ model of the perturbing potential, we compared our Figure~\ref{fig:Fig599} with figure 2 of \cite{Abad}, obtaining the same result. Now, using the extended $J_{2}+J_{3}+\cdots+J_{7}$ model, our Figures~\ref{fig:Fig601},~\ref{fig:Fig602}, and~\ref{fig:Fig603} agree with figures 15(b) and 17(a) of \cite{Sirwah}. This confirms that our approach is suitable for testing with a slightly more robust model, that is, considering the zonal terms up to the $J_{18}$ term plus the third-body perturbation. Note that in both structures, we obtain frozen orbits in regions with inclinations close to 90 degrees, except in Figure~\ref{fig:Fig599}, where the eccentricity variation amplitude is very large, leading to an impact of the satellite with the lunar surface. This behavior occurs because low-altitude orbits are strongly affected by the lunar gravitational field.

Finally, considering our most complete model, referred to as Model $18 \times 0$, we present Figures~\ref{fig:Fig604} and~\ref{fig:Fig605}. New orbital solutions appear; however, they exhibit high eccentricities, making them impractical for medium- and low-altitude missions. Figure~\ref{fig:Fig605} shows that the frozen-orbit region near $90$ degrees inclination persists, but now corresponds to an initial value of $g = 270$ degrees, whereas in the previous figures these orbits were obtained for $g = 90$ degrees. To conclude this approach, we consider low-altitude orbits around the Moon. Figure~\ref{fig:Fig606} ($a = 1838$ km) shows a polar ($i=90$ degrees or 1.57 radians) and frozen orbit with $e = 0.049$ {\color{black}(this value will be discussed later)}, while Figure~\ref{fig:Fig608} ($a = 1938$ km) presents a polar frozen orbit with $e = 0.0241$. Additionally, Figure~\ref{fig:Fig608} reveals the existence of a frozen orbit with inclination $78.2087$ degrees and eccentricity $0.014$, with $g = 90$ degrees, which may be of interest for scientific missions around the Moon.
   
In this approach, considering Model $18 \times 0$ within Lagrange's planetary equations and imposing the condition $dg/dt = 0$, instead of deriving an analytical expression of the form $e = f(i, g)$, which becomes excessively complex, we performed a numerical simulation (more precisely, a numerical sweep to identify equilibrium points) and subsequently plotted the numerically obtained results. This procedure generates a numerical representation of the equilibrium points, expressed as $e = f(i, a, g)$, under the condition $dg/dt = 0$ while fixing the values of $a$ and $g$. According to \cite{Kikuchi}, frozen-orbit solutions under irregular gravity perturbations and solar radiation pressure can be obtained by numerically solving the equation $dg/dt = 0$. Since the determination of these solutions does not require full numerical propagation, the proposed orbit design method can be classified as semi-analytical. Using this methodology, we constructed Figures~\ref{fig:Fig506} and~\ref{fig:Fig508}.  


We analyze the frozen orbit conditions considering the selenopotential model, namely, Model 18 $\times$ 0 ($J_{2}+\cdot \cdot \cdot +J_{18}$).	According to Figure ~\ref{fig:Fig605}, we considered only the value of $g=270$ degrees, since only this value contributed to obtaining a frozen polar orbit. Figures ~\ref{fig:Fig506} and ~\ref{fig:Fig508} show that for both values of semi-major axis, frozen orbits were obtained. But in the case of $a=1838$ km, the eccentricity value of the frozen polar orbit is very high, $e=0.049$ (Note that this is consistent with the Figure ~\ref{fig:Fig605}), and as we will see in the next section, this orbit does not survive for very long. For the case of the semi-major axis of 1938 km, we find a frozen polar orbit, $e=0.024$, which, as we will see in the next section, survives for quite some time. Note that these figures also show other regions of frozen orbits and regions of orbits with large variations in eccentricity, especially around the 60 and 120 degrees inclination. But at the moment, we are interested in the frozen polar orbits, which are generally of interest for lunar science missions. The inclination region near 60 and 120 degrees will also be discussed in the resonance section.
\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig506.eps}
		\caption{Frozen conditions for an orbiter: $a=1838$ km, $e\in \{0.001,0.3\}$, $g=270$ deg, $h=90$ deg. Model considering only the zonal harmonics $J_{2}+\cdot\cdot\cdot+J_{18}$. Frozen polar orbit: $e=0.049$.}
		\label{fig:Fig506}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig508.eps}
		\caption{Frozen conditions for an orbiter: $a=1938$ km, $e\in \{0.001,0.3\}$, $g=270$ deg, $h=90$ deg. Model considering only the zonal harmonics $J_{2}+\cdot\cdot\cdot+J_{18}$. Frozen polar orbit: $e=0.024$.}
		\label{fig:Fig508}
	\end{minipage}
\end{figure}



\subsection{Numerical simulation: frozen orbits}

In this section, we select the frozen initial conditions obtained in the analytical study. These initial conditions are propagated using a numerical integrator that incorporates a potential function taking into account, in addition to the zonal terms, the sectoral and tesseral terms of order and degree 3, given by:
\begin{equation}\label{potential}
	\begin{array}{l}
		\displaystyle \langle R \rangle = \langle \sum _{n=2}^{18} R_{J_{n}}\rangle + \langle R_{C_{22}} \rangle + \langle R_{S_{22}} \rangle + \langle \sum _{m=1}^{3} R_{C_{3m}}  \rangle + \langle \sum _{m=1}^{3} R_{S_{3m}}  \rangle.
	\end{array}
\end{equation}
The disturbing potential given in Eq. (\ref{potential}) is replaced in the Lagrange planetary equations and numerically integrated to investigate the behavior of the orbital elements of artificial satellites and/or space debris. This system of equations is numerically integrated using Maple software. For the numerical integration of the dynamical system, we use Maple dsolve routine with the options numeric and method = rkf45, which represents an integration using the Runge-Kutta-Fehlberg method. In this section, we analyze the behavior of low-altitude polar orbits around the Moon. So, numerically integrating equation (\ref{potential}), via Lagrange's planetary equations, we obtain Figures  ~\ref{fig:Fig505} and ~\ref{fig:Fig507}.

Figure~\ref{fig:Fig505} shows the diagram of eccentricity versus the argument of perilune for the frozen orbit, obtained from Figures~\ref{fig:Fig506} and~\ref{fig:Fig606}, but now through numerical simulation. Different eccentricity values are considered, including the obtained value $e = 0.049$, using Model $18 \times 0$. Note that the frozen orbit is highly eccentric for the considered semi-major axis ($1838$ km). Therefore, even a small perturbation is sufficient for an artificial satellite in this orbit to impact the lunar surface very quickly. In this case, eccentricity values other than that of the frozen orbit ($e = 0.049$) are also considered. These additional eccentricity values cause the orbits to librate around the equilibrium point with a certain amplitude of variation.

In fact, when we perform the numerical simulation considering our most complete model ( $\langle {J_{2}+\cdot\cdot\cdot+J_{18}\rangle}+\langle {C_{22}} \rangle + \langle {S_{22}} \rangle + \langle \sum _{m=1}^{3} {C_{3m}}  \rangle + \langle \sum _{m=1}^{3} {S_{3m}}  \rangle$), Model $18 \times$ 3, we obtain Figure ~\ref{fig:Fig507}. Note that the frozen orbit begins to have a small amplitude of variation, which we call a quasi-frozen orbit. However, because it is very eccentric, even with a small variation in the eccentricity and perilune argument, the satellite in this orbit impacts in a little over 250 days.



\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig505.eps}
		\caption{Initial conditions: $a=1838$ km, $i=90.295$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.049$. Model considering only the zonal harmonics, Model 18 $\times$ 0.}
		\label{fig:Fig505}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig507.eps}
		\caption{Initial conditions: $a=1838$ km, $i=90.295$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.049$. Model considered, Model 18 $\times$ 3.}
		\label{fig:Fig507}
	\end{minipage}
\end{figure}

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig610.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90.2981$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.024$. Model considering only the zonal harmonics, Model 18 $\times$ 0.}
		\label{fig:Fig610}
	\end{minipage}%
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Fig611.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90.2981$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.024$. Model considered, Model 18 $\times$ 3.}
		\label{fig:Fig611}
	\end{minipage}
\end{figure}





Now, considering the frozen orbit obtained in the analytical study with a semi-major axis of $1938$ km plotted in Figures~\ref{fig:Fig608} and~\ref{fig:Fig508}, we performed numerical integrations using Maple, according to the models described in the figure legend. In this case, we also consider our most complete model. As in the previous case, we reproduce the frozen orbit ($e = 0.024$) and include additional eccentricity values. As Figures~\ref{fig:Fig610} and \ref{fig:Fig611} show, the orbital elements of the frozen orbit (see Figure~\ref{fig:Fig610}) have their amplitude increased when our more complete model is considered (see Figure~\ref{fig:Fig611}). Thus, we obtain a quasi-frozen orbit, that is, an orbit that presents less variation in the orbital elements, which can be evaluated in an even more precise model, such as, for example, with the inclusion of orbital perturbation from other planets, relativistic effects, among others. With the other eccentricity values considered, the orbits are strongly perturbed.
    
  
  

Finally, Figures ~\ref{fig:grafico1} and ~\ref{fig:grafico3} show the quasi-frozen orbit obtained in Figure ~\ref{fig:Fig611}, where $e=0.0241$ and $a=1938$ km. Figure ~\ref{fig:grafico3} shows a comparison of our model with a model of order and degree $50 \times 50$ (using the GMAT software). Considering our Model 18 $\times$ 3, we show that it represents the behavior of the orbit with good precision, although it does not faithfully represent Model 50 $\times$ 50, it approximates it reasonably well for a first analysis (see Figure ~\ref{fig:grafico3}). Thus, the characteristic of the orbit remains, that is, the quasi-frozen orbit remains even with the inclusion of a more complete model, considering the high order of spherical harmonics. Note that the orbit found presents small variations in eccentricity and perilune argument (see Figure ~\ref{fig:grafico1}), which we call a quasi-frozen orbit.

Therefore, the initial conditions for a quasi-frozen polar orbit, obtained in this work, for a low-altitude satellite around the Moon are: $a=1938$ km, $e=0.0241$, $i=90$ degrees, $g=270$ degrees and $h=90$ degrees. We suggest that these data be evaluated in more complex models to verify the possibility of maintaining the quasi-frozen polar orbit.

\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{grafico1.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.0241$. GMAT Model 50 $\times$ 50}
		\label{fig:grafico1}
	\end{minipage}
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{grafico3.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.0241$. GMAT Model 50 $\times$ 50, and Semi-analytical Model 18 $\times$ 3.}
		\label{fig:grafico3}
	\end{minipage}%
\end{figure}

To find the frozen orbit in the analytical model, we set $g=270$ degrees. In the numerical simulations, we used a random value for the longitude of the ascending node ($h=90$ degrees). Now, we show a dynamic map of the argument of the perilune versus the longitude of the ascending node to verify the contribution of these orbital elements to the increase in eccentricity. We used a 5-degree step for $g$ and $h$. Thus, Figure~\ref{fig:grafico15} shows that, in fact, the value of $g$ close to 270 degrees contributes to a smaller amplitude of eccentricity variation, as we can see in the color scale in Figure~\ref{fig:grafico15}. Note that for other values of $g$, there is a large variation in eccentricity. As we can see in Figure~\ref{fig:grafico15}, there is practically no contribution from the longitude of the ascending node to the increase in eccentricity.



\begin{figure}[!h]
	\centering
	\resizebox{0.7\columnwidth}{!}{\includegraphics{grafico15.eps} }
	\caption{Initial conditions: $a=1938$ km, $i=90$ deg. Frozen polar orbit: $e=0.0241$. $g=0..360$ deg ($\Delta g=5$ deg),  $h=0..360$ deg ($\Delta h=5$ deg). Integration time: 3 years. Model 18 $\times$ 3.}
	\label{fig:grafico15}       % Give a unique label
\end{figure}


\begin{figure}[!h]
	\centering
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Grafico40.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.0241$. GMAT Model 50 $\times$ 50 + Earth, and Semi-analytical Model 18 $\times$ 3 + Earth.}
		\label{fig:Grafico40}
	\end{minipage}
	\hfill
	\begin{minipage}[t]{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth,height=5.5cm]{Grafico41.eps}
		\caption{Initial conditions: $a=1938$ km, $i=90$ deg, $g=270$ deg, $h=90$ deg. Frozen polar orbit: $e=0.0241$. GMAT Model 50 $\times$ 50 + Earth, and Semi-analytical Model 18 $\times$ 3 + Earth.}
		\label{fig:Grafico41}
	\end{minipage}
\end{figure}

To conclude this section, in addition to the orbital perturbations given in Eq.~(\ref{potential}), we include the third-body perturbation due to the Earth in order to make the model more accurate. In the semi-analytical model, the disturbing potential is given by \cite{Carvalho3} (see Appendix B), where the Earth's orbit is assumed to be elliptical and inclined, adopting a single averaged formulation. For comparison purposes, we use the Earth's gravity model implemented in the GMAT software. The results show that the inclusion of the third-body perturbation only slightly modifies the eccentricity variation as it is illustrated in Figures~\ref{fig:grafico1},~\ref{fig:grafico3},~\ref{fig:Grafico40}, and~\ref{fig:Grafico41}. Therefore, the obtained orbit remains quasi-frozen in both models.



\subsection{Analytical and numerical model: resonant orbits}
An analytical approach is conducted to study the six resonance cases presented in the previous section (see Eq.~\ref{eq:resonant_equations}). The goal is to determine initial values of the semi-major axis and inclination that contribute to amplifying the rate of change of orbital eccentricity and inclination, which may accelerate the reentry of space debris into the lunar surface. On the other hand, we also seek values of the semi-major axis and inclination that stabilize eccentricity and inclination, thereby helping to delay the reentry of a spacecraft.

For this purpose, we consider a single solar radiation pressure (SRP) resonance term, $\langle R_{SRP_i} \rangle$, and apply Lagrange’s planetary equations to compute the variations of the orbital elements. Without loss of generality, we analyze the resonance $SRP_1$. Under this assumption, the eccentricity evolution, when only this resonance is considered, is given by

\begin{equation}\label{eqLagwe}
	\frac{de}{dt} = \frac{8\sqrt{1-e^2}}{27\,n\,a}\,(1+\cos i)(1+\cos \epsilon)\frac{\partial \cos \psi_1}{\partial g},
\end{equation}

where the resonance term $\dot{\psi}_{1}$ introduces a singularity in Eq.~(\ref{eqLagwe}) during integration. Similarly, the following expressions appear when analyzing the other resonances:

\begin{equation}\label{resonance}
	\begin{array}{rlclrll}
		\dot{\psi}_{1}&=&\dot{h}+\dot{g}-\dot{\lambda_{\odot}},\qquad
		\dot{\psi}_{2}&=&\dot{h}-\dot{g}-\dot{\lambda_{\odot}},\qquad
		\dot{\psi}_{3}&=&\dot{g}-\dot{\lambda_{\odot}},\\
		\dot{\psi}_{4}&=&\dot{g}+\dot{\lambda_{\odot}},\quad\qquad
		\dot{\psi}_{5}&=&\dot{h}+\dot{g}+\dot{\lambda_{\odot}},\qquad
		\dot{\psi}_{6}&=&\dot{h}-\dot{g}+\dot{\lambda_{\odot}}.
	\end{array}
\end{equation}

Here, $\dot{h}$, $\dot{g}$, and $\dot{\lambda_{\odot}}$ denote time derivatives. The quantities $\dot{h}$ and $\dot{g}$ represent the precession rates of the right ascension of the ascending node and the argument of periselene, respectively, while $\dot{\lambda_{\odot}} = n_{\odot}$ since $\lambda_{\odot} = n_{\odot} t + \lambda_{0}$, where $\lambda_{0}$ is a constant defined at the initial epoch.

At this stage, in order to obtain analytical insight into each resonance through Lagrange’s planetary equations, we assume that the disturbing potential is given by

\[
R = \langle R_{J_{2}} \rangle + \langle R_{SRP} \rangle .
\]

This assumption is adopted because $J_{2}$ represents the dominant term of the lunar gravitational field, while SRP is the main non-gravitational perturbation responsible for generating the resonances under investigation. Therefore, considering only these contributions allows the derivation of simplified analytical relations that capture the primary resonant dynamics.

However, the lunar gravitational field differs significantly from the terrestrial case. In particular, the lunar coefficient $J_{7}$ differs from $J_{2}$ by only about one order of magnitude, whereas for the Earth this difference is much larger. 

Consequently, a more realistic lunar model should include higher-degree harmonics. When these additional terms are incorporated, the resonance condition leads to an expression for the semi-major axis that depends not only on eccentricity and inclination, but also on the argument of periselene, resulting in a relation of the form $a = f(e,i,g)$. This dependence implies that the semi-major axis varies with the argument of periselene, making the analytical treatment of the resonance significantly more complex. Moreover, different values of $g$ generate distinct families of resonant orbits. For this reason, we restrict the analytical development to the simplified model described above.



It is possible to obtain the time variation of $g$ and $h$ for each resonance given in Eq.~(\ref{eq:resonant_equations}) as follows:

\begin{equation}\label{28}
\frac{dg}{dt}=\frac{dg_{SRP_{i}}}{dt}+\frac{dg_{J_{2}}}{dt},
\end{equation}

\begin{equation}\label{29}
\frac{dh}{dt}=\frac{dh_{SRP_{i}}}{dt}+\frac{dh_{J_{2}}}{dt},
\end{equation}
where $i=1,\ldots,6$. By imposing the resonance conditions given in Eq.~(\ref{resonance}) equal to zero, we obtain an expression to determine the resonant semi-major axis associated with each SRP term, which becomes a function of inclination and eccentricity, i.e., $a(i,e)$. This relation allows us to identify the location of each resonance for selected values of eccentricity and area-to-mass ratio.

The following constant parameters were adopted: $A/m=0.012$ m$^{2}$/kg and $\lambda_{\odot}=90.086$ degrees, corresponding to the initial epoch of June 21 at 06:43:12, as taken from \citet{Alessi}. By fixing the eccentricity and varying the inclination, the resonant semi-major axis for each resonance can be obtained, as shown in Figures~\ref{fig:Fig64} and~\ref{fig:Fig65}. Figure~\ref{fig:Fig64} presents the resonant locations without considering the physical radius of the Moon, thus including regions that do not correspond to realistic orbital motion because they represent trajectories below the lunar surface. In contrast, Figure~\ref{fig:Fig65} displays only the resonant regions associated with physically valid semi-major axis values, corresponding to low-altitude lunar orbits located above the lunar radius.

\begin{figure}[!h]
\centering
\resizebox{0.5\columnwidth}{!}{\includegraphics{Fig64.eps}}
\caption{$(a,i)$ plot of all types of resonant orbits for dynamical model $J_{2}+ SRP$. Initial conditions: $e=0.01$, $g=270^\circ$, $h=90^\circ$, $A/m=0.012$ m$^{2}/$kg.}
\label{fig:Fig64}
\end{figure}

\begin{figure}[!h]
\centering
\resizebox{0.5\columnwidth}{!}{\includegraphics{Fig65.eps}}
\caption{$(a,i)$ plot of all types of resonant orbits for dynamical model $J_{2}+ SRP$. Initial conditions: $e=0.01$, $g=270^\circ$, $h=90^\circ$, $A/m=0.012$ m$^{2}/$kg. Physically valid semi-major axis values, corresponding to low-altitude lunar orbits.}
\label{fig:Fig65}
\end{figure}

To analyze the combined effect of all resonances shown in Figure~\ref{fig:Fig65}, numerical simulations were performed considering only SRP and the $J_{2}$ contribution in the disturbing potential (see Figure~\ref{fig:grafico10}). The dynamical map presented in Figure~\ref{fig:grafico10} illustrates the amplitude of eccentricity variation, defined as $Amp(e)=e_{max}-e_{min}$. The results indicate that the resonance effect produces only modest eccentricity variations. A small region corresponding to very low inclinations and altitudes exhibits an increase in eccentricity for inclinations below approximately $20^\circ$. This behavior differs from the Earth case, where the coupling between the $J_{2}$ term and SRP can strongly affect eccentricity for specific inclination values (see \citet{Alessi}).

\begin{figure}[!h]
\centering
\begin{minipage}[t]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth,height=5.5cm]{grafico10.eps}
\caption{Dynamic map of $Amp(e)=e_{max}-e_{min}$ over 5 yr. $(a,i)$ plot of all types of resonant orbits for dynamical model $J_{2}+ \langle U_{SRP}\rangle$. Initial conditions: $a=1800..2700$ km ($\Delta a=10$ km), $e=0.01$, $i=1..175^\circ$ ($\Delta i=1^\circ$), $g=270^\circ$, $h=90^\circ$, $A/m=0.012$ m$^{2}/$kg.}
\label{fig:grafico10}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth,height=5.5cm]{grafico14.eps}
\caption{$(a,i)$ plot of all types of resonant orbits for dynamical model Model $18 \times 0 + \langle U_{SRP}\rangle$. Same initial conditions as Figure~\ref{fig:grafico10}.}
\label{fig:grafico14}
\end{minipage}
\end{figure}

To better understand the coupling between the lunar gravitational field and SRP, additional zonal harmonics were included, extending the model up to the $J_{18}$ term. Figure~\ref{fig:grafico14} presents the resulting dynamical behavior. The inclusion of higher-degree harmonics significantly modifies the dynamical structure, producing larger eccentricity amplitudes and generating additional resonance regions. This effect is particularly evident near inclinations of approximately $60$ degrees and $120$ degrees, which correspond to the critical inclination region for prograde and retrograde orbits. As shown in Figure~\ref{fig:grafico14}, the semi-major axis interval between 2000 km and 2300 km presents a strong eccentricity growth for resonant inclinations close to these critical values. The polar region, on the other hand, is mainly affected for low-altitude orbits up to approximately 1900 km.

Considering the same initial conditions used in Figure~\ref{fig:grafico14}, we further expanded the dynamical model by including sectoral and tesseral harmonics, as well as the gravitational perturbation of the Earth. This configuration corresponds to Model $18 \times 3 + \langle U_{SRP} \rangle + \langle U_{Earth} \rangle$, where harmonics up to order and degree three are considered. The results are presented in Figure~\ref{fig:grafico11}. The inclusion of these additional harmonic coefficients produces substantial changes in the eccentricity behavior, as can be observed by comparing Figures~\ref{fig:grafico10},~\ref{fig:grafico14}, and~\ref{fig:grafico11}. Due to the irregular nature of the lunar gravitational field, the orbital dynamics are strongly modified when higher-order terms are included.

\begin{figure}[!h]
\centering
\begin{minipage}[t]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth,height=5.5cm]{grafico11.eps}
\caption{Dynamic map of $Amp(e)=e_{max}-e_{min}$ over 5 yr. $(a,i)$ plot of all types of resonant orbits for dynamical model Model $18 \times 3 + \langle U_{SRP}\rangle + \langle U_{Earth}\rangle$. Initial conditions: $a=1800..2700$ km ($\Delta a=10$ km), $e=0.01$, $i=1..175^\circ$ ($\Delta i=1^\circ$), $g=270^\circ$, $h=90^\circ$, $A/m=0.012$ m$^{2}/$kg.}
\label{fig:grafico11}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth,height=5.5cm]{grafico12.eps}
\caption{Dynamic map of the satellite lifetime over 5 years. Same initial conditions as Figure~\ref{fig:grafico11}.}
\label{fig:grafico12}
\end{minipage}
\end{figure}

In regions close to the critical inclinations (both prograde and retrograde), the resonance mechanism substantially increases eccentricity growth. This effect may be useful for end-of-life disposal strategies, facilitating the removal of spacecraft from lunar orbit, but may also represent a risk for operational satellites located within these resonant regions. It is important to note that a typical spacecraft with $A/m=0.012$ m$^{2}$/kg was considered, implying that SRP produces only a moderate dynamical effect.

Figure~\ref{fig:grafico12} shows the satellite lifetime represented through a color scale. The results indicate that low-altitude orbits, including polar trajectories, tend to collide with the lunar surface within relatively short timescales for several inclination values. Orbits near the critical inclination regions exhibit amplified eccentricity growth, affecting a broader range of semi-major axis values. As expected, higher-altitude debris presents longer collision times. From Figure~\ref{fig:grafico12}, it can be observed that the polar orbit with a semi-major axis of 1938 km remains stable and does not collide with the lunar surface within the analyzed time span.

The simulations also included the third-body perturbation due to the Earth. The results show behavior similar to that obtained with Model $18 \times 3 + \langle U_{SRP} \rangle$, indicating that the overall structure of the dynamical map remains essentially unchanged. The Earth’s perturbation produces only small increases in eccentricity, slightly expanding the affected semi-major axis regions. Nevertheless, the region near the semi-major axis of 1938 km remains practically unaffected.

\section{Conclusions}

In this work, we investigate the dynamics of artificial satellites orbiting the Moon considering the main disturbing forces. In particular, we analyze the orbits under the effect of the lunar gravitational potential taking into account a simplified model, considering up to the order eighteen for the zonal harmonics $J_{2}+\cdot \cdot \cdot +J_{18}$, and the sectoral and tesseral terms of order and degree three. Considering our simplified model (Model 18 $\times$ 3 $+ \langle U_{SRP}\rangle+\langle U_{Earth}\rangle$), we investigate long and short duration orbits around the Moon for different types of missions, namely, low-altitude artificial satellites and space debris disposal. 
	
With the analytical model, we analyze the conditions of frozen orbits when the perilune argument ($g$) is 270 and 90 degrees. We show the results in graphs of eccentricity versus inclination, where we obtain initial conditions of frozen polar orbits for $g=270$ degrees. Numerically integrating the equations of motion using the initial conditions obtained in the analytical model, we identify quasi-frozen polar orbits. After analyzing the results, we find a frozen polar orbit of low altitude given by the initial conditions: $a=1938$ km, $e=0.0241$, $i=90$ degrees, $g=270$ degrees, and $h=90$ degrees. To validate the results, we used the GMAT software considering the selenopotential model of order and degree 50 $\times$ 50. Considering our Model 18 $\times$ 3, we show that it represents the behavior of the orbit with good accuracy when compared to Model 50 $\times$ 50 (GMAT). We show a dynamic map of the argument of the perilune versus longitude of the ascending node to verify the contribution of these orbital elements to increased eccentricity. In fact, the value of $g$ close to 270 degrees contributes to a lower range of eccentricity variation. We observed that for other values of $g$, there is a large variation in eccentricity. We note that there is practically no contribution of the longitude of the ascending node to increased eccentricity.
	
We find that, in general, the resonance effect due to the coupling of the $J_{2}$ term and the solar radiation pressure (SRP) does not contribute to amplify the growth of the eccentricity, except for some low inclination values. Due to the characteristics of the lunar gravitational field, we included the zonal harmonics up to the $J_{18}$ term. With the inclusion of these coefficients, we obtain a large variation in the behavior of the eccentricity. Due to the characteristics of the selenopotential, the dynamics change completely when more terms are added. To better understand the resonance effect due to the coupling of the lunar gravitational field and the SRP, we include zonal harmonics up to order eighteen and sectoral and tesseral terms of order and degree three. Using the dynamic map of the semi-major axis versus inclination, we show the regions where eccentricity is strongly affected, where more resonance zones emerge. These regions of large eccentricity variation can aid in the removal of space debris from lunar orbit. In particular, we show a region of the semi-major axis between 2000 km and 2300 km that exhibits a significant increase in eccentricity for resonant inclinations close to 60 and 120 degrees (where a clear resonance phenomenon emerges), corresponding to the critical inclination region in prograde and retrograde cases. We also show the lifetime of the satellite and/or space debris when placed in the resonance regions. Using the dynamic map, we can also identify the regions where an artificial satellite in lunar orbit should avoid due to resonance effects. Another important point is when the effect of the third-body perturbation (Earth) is taken into account. The result shows that the characteristic dynamic map remains, with only small amplifications of the eccentricity being found. 



%%%%%%%%%%%%%%%%%%%%%%%%%%


\backmatter

%\bmhead{Supplementary information}
%
%If your article has accompanying supplementary file/s please state so here. 
%
%Authors reporting data from electrophoretic gels and blots should supply the full unprocessed scans for key as part of their Supplementary information. This may be requested by the editorial team/s if it is missing.
%
%Please refer to Journal-level guidance for any specific requirements.

%\bmhead{Acknowledgements}
%
%Acknowledgements are not compulsory. Where included they should be brief. Grant or contribution numbers may be acknowledged.
%
%Please refer to Journal-level guidance for any specific requirements.

\section*{Declarations}

%Some journals require declarations to be submitted in a standardised format. Please check the Instructions for Authors of the journal to which you are submitting to see if you need to complete this section. If yes, your manuscript must contain the following sections under the heading `Declarations':

\begin{itemize}
\item Funding. Sponsored by CNPq - Brazil (National Council for Scientific and Technological Development) for contract 408221/2025-9 and FAPESP: proc. 2022/11783-5. This work has been partially supported by Grant PID2024-156002NB-I00 by MICIU$/$AEI$/$10.13039$/$501100011033$/$FEDER, UE, and by the Arag\'on Government and European Social Fund (E24$\_$23R).
\item AI tools were used to improve the English language and clarity of the manuscript.
\item Conflict of interest/Competing interests. Not applicable. 
\item Ethics approval and consent to participate. Not applicable. 
\item Consent for publication. Not applicable. 
\item Data availability. The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. 
\item Materials availability. Not applicable. 
\item Code availability. The different codes used in this paper can be reproduced since all the formulations are included. 
\item Author contribution. All authors have contributed equally in the present manuscript. All authors have read and approved the final manuscript. 
\end{itemize}


%\noindent
%If any of the sections are not relevant to your manuscript, please include the heading and write `Not applicable' for that section. 

%%===================================================%%
%% For presentation purpose, we have included        %%
%% \bigskip command. Please ignore this.             %%
%%%===================================================%%
%\bigskip
%\begin{flushleft}%
%Editorial Policies for:
%
%\bigskip\noindent
%Springer journals and proceedings: \url{https://www.springer.com/gp/editorial-policies}
%
%\bigskip\noindent
%Nature Portfolio journals: \url{https://www.nature.com/nature-research/editorial-policies}
%
%\bigskip\noindent
%\textit{Scientific Reports}: \url{https://www.nature.com/srep/journal-policies/editorial-policies}
%
%\bigskip\noindent
%BMC journals: \url{https://www.biomedcentral.com/getpublished/editorial-policies}
%\end{flushleft}

\begin{appendices}
	
	%\section{Section title of first appendix}\label{secA1}
	
	
	\section{Averaged disturbing potential terms due to the non-uniform distribution of mass of the central body} \label{secA1}
	
	In \citet{Coffey}, a model of Earth's gravitational field is presented, expanded up to the $J_{9}$ term. However, the authors write the equations normalized by $J_{2}^{2}$ and use some parameters that make comparison difficult in an initial analysis. Therefore, we did not compare the models.
	
	\begin{equation*}\label{6}
		\begin{array}{l}
			\langle R_{J2} \rangle=-\frac{1}{4}\epsilon{\frac {\sqrt {5}{n}^{2} \left( 3\,{s}^{2}-2 \right) }
				{ \left( -{e}^{2}+1 \right) ^{3/2}}},		
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{7}
		\begin{array}{l}
			\langle R_{J3} \rangle=-{\frac {15}{8}{\epsilon_{1}}}{\frac {\sqrt {7}{n}^{2} e s \left( {s}^
					{2}-4/5 \right) \sin \left( g \right) }{ \left( -{e}^{2}+1 \right) ^
					{5/2}a}},
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{8}
		\begin{array}{l}
			\langle R_{J4} \rangle={\frac {105}{64}}{\epsilon_{2}}\frac{\sqrt{9}{n}^{2}}{ \left( -{e}^{2}+1 \right) ^{7/2}{a}^{2}} \left( {e}^{2}{s}^{2} \left( {s}^{2}-6/7 \right) \cos \left( 2\,g \right) -3/2\, \left( {s}^{4}-{\frac {
					8}{7}}\,{s}^{2}+{\frac {8}{35}} \right)  \left( {e}^{2}+2/3 \right)
			\right),		
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{18}
		\begin{array}{l}
			\langle R_{J5} \rangle=-{\frac {945}{128}}{\epsilon_{3}}\frac{\sqrt {11}{n}^{2} e s}{ \left( -{e}
				^{2}+1\right)^{9/2}{a}^{3}} \left( -1/6\, \left( {s}^{2}-{
				\frac {8}{9}} \right) {e}^{2}{s}^{2}\sin \left( 3\,g \right) + \left(
			{e}^{2}+4/3 \right)  \left( {s}^{4}-4/3\,{s}^{2}+{\frac {8}{21}}
			\right) \sin \left( g \right)  \right),
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{19}
		\begin{array}{ll}
			\langle R_{J6} \rangle & ={\frac {17325}{2048}}{\epsilon_{4}}\frac{\sqrt {13}{n}^{2}}{( -{e}^{2}+1)^{11/2}{a}^{4}}(( {s}^{4}-{\frac
				{16}{11}}{s}^{2}+{\frac {16}{33}}) {e}^{2}{s}^{2}( {e}^
			{2}+2) \cos ( 2g) -1/10{e}^{4}{s}^{4} ( {s
			}^{2}-{\frac {10}{11}} ) \cos( 4g) \\[1.5ex]
			\qquad & - ( {e}^{
				4}+8/3{e}^{2}+{\frac {8}{15}})( {s}^{6}-{\frac {18}{11
			}}{s}^{4}+{\frac {8}{11}}{s}^{2}-{\frac {16}{231}})), \\
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{20}
		\begin{array}{ll}
			\langle R_{J7} \rangle & =-{\frac {9009}{16384}}{\epsilon_{5}}\frac{\sqrt {15}{n}^{2} e s}{( -{e}^{2}+1)^{13/2}{a}^{5}}( -15( {
				s}^{4}-{\frac {20}{13}}{s}^{2}+{\frac {80}{143}}) {s}^{2}{e}^
			{2}( {e}^{2}+8/3) \sin( 3g) +{s}^{4}( {s}^{2}-{\frac {12}{13}}) {e}^{4}\sin( 5g) \\[1.5ex]
			\qquad & + 50( {s}^{6}-{\frac {24}{13}}{s}^{4}+{\frac {144}{
					143}}{s}^{2}-{\frac {64}{429}})( {e}^{4}+4{e}^{2}+8/
			5) \sin( g ) ),
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{21}
		\begin{array}{ll}
			\langle R_{J8} \rangle & = {\frac {45045}{131072}}{\epsilon_{6}}\frac{\sqrt {17}{n}^{2}}{( -{e}^{2}+1)^{15/2}{a}^{6}}( 105( {e}^{4}+16/3{e}^{2}+{
				\frac {16}{5}})( {s}^{6}-2{s}^{4}+{\frac {16}{13}}{s
			}^{2}-{\frac {32}{143}}) {e}^{2}{s}^{2}\cos( 2g) \\[1.5ex]
			\qquad & - 21{e}^{4}( {e}^{2}+10/3) {s}^{4}( {s}^{
				4}-8/5{s}^{2}+{\frac {8}{13}}) \cos ( 4g ) +{e}^
			{6}{s}^{6}( {s}^{2}-{\frac {14}{15}}) \cos( 6g) \\[1.5ex]
			\qquad & - {\frac {175}{2}}( {s}^{8}-{\frac {32}{15}}{s}^{6}+
			{\frac {96}{65}}{s}^{4}-{\frac {256}{715}}{s}^{2}+{\frac {128}{
					6435}})( {e}^{6}+6{e}^{4}+{\frac {24}{5}}{e}^{2}+{
				\frac {16}{35}})),
		\end{array}
	\end{equation*}
	
	\begin{equation*}\label{22}
		\begin{array}{ll}
			\langle R_{J9} \rangle & = -{\frac {765765}{131072}}{\epsilon_{7}} \frac{\sqrt {19}{n}^{2}e s}{( -{e}^{2}+1)
				^{17/2}{a}^{7}}( -7\,{e}^{2}
			( {s}^{6}-{\frac {36}{17}}\,{s}^{4}+{\frac {24}{17}}\,{s}^{2}-{
				\frac {64}{221}}) {s}^{2}( {e}^{4}+{\frac {20}{3}}\,{e}^{
				2}+16/3 ) \sin ( 3g ) \\[1.5ex]
			\qquad & + {e}^{4}{s}^{4}( {e}^{2
			}+4 )( {s}^{4}-{\frac {28}{17}}\,{s}^{2}+{\frac {56}{85}}) \sin ( 5\,g ) -1/28\,{e}^{6} ( {s}^{2}-{
				\frac {16}{17}}) {s}^{6}\sin ( 7\,g ) \\[1.5ex]
			\qquad & + {\frac {35}{2
			}}\,\sin( g)( {\frac {64}{35}}+{e}^{6}+8\,{e}^{4}+
			{\frac {48}{5}}\,{e}^{2})( {s}^{8}-{\frac {40}{17}}\,{s}
			^{6}+{\frac {32}{17}}\,{s}^{4}-{\frac {128}{221}}\,{s}^{2}+{\frac {128
				}{2431}})),
		\end{array}
	\end{equation*}
	where $\epsilon=J_{2}R_{M}^{2}$, $\epsilon_{1}=J_{3}R_{M}^{3}$, $\epsilon_{2}=J_{4}R_{M}^{4}$, $\epsilon_{3}=J_{5}R_{M}^{5}$, $\epsilon_{4}=J_{6}R_{M}^{6}$, $\epsilon_{5}=J_{7}R_{M}^{7}$, $\epsilon_{6}=J_{8}R_{M}^{8}$, $\epsilon_{7}=J_{9}R_{M}^{9}$, $c=\cos(i)$ and $s=\sin(i)$.
	
	
	\begin{equation*}\label{23}
		\begin{array}{l}
			\langle R_{C22,S22} \rangle   =\frac{3}{2}\sqrt {\frac{5}{12}} \frac {{n}^{2}{R}^{2}_{M} {s}^{2} }{ \left(1 -{e}^{2}\right)^{3/2}}
			\left( {C_{22}}\,\cos \left( 2\,{\nu}\,t-2\,h \right) -{S_{22}}
			\,\sin \left( 2\,{\nu}\,t-2\,h \right)  \right),
		\end{array}
	\end{equation*}
	
	
	\begin{equation*}\label{24}
		\begin{array}{ll}
			\langle R_{C31,S31} \rangle  &= -{\frac {15}{8}}\sqrt {\frac{7}{6}} \frac{ e {n}^{2}{R}^{3}_{M}}{( 1-{e}^{2}
				)^{5/2}{a}}(( {C_{31}}( {s
			}^{2}-4/5) \cos ( g) +3c( {s}^{2}-{\frac {4}
				{15}}) \sin ( g) {S_{31}} ) \cos( {\nu}\,t-h)\\[1.5ex]
			\qquad & + 3 \sin( {\nu}\,t-h)(
			( -1/3\,{s}^{2}+{\frac {4}{15}}) {S_{31}}\,\cos ( g
			) + ( {s}^{2}-{\frac {4}{15}}) {C_{31}}\,\sin( g) c)),
		\end{array}
	\end{equation*}
	
	
	\begin{equation*}\label{25}
		\begin{array}{ll}
			\langle R_{C32,S32} \rangle  & =\frac{15}{2} \sqrt {\frac{7}{120}} \frac {e s{n}^{2}{R}^{3}_{M}}{( 1-{e}^
				{2} )^{5/2}a}((( -3/2{C_{32}}{c}^{2}+1/2\,{
				C_{32}}) \sin( g) +\cos( g) { S_{32}
			}\,c) \cos ( 2\,{\nu}\,t-2\,h ) \\[1.5ex]
			\qquad & + \sin ( 2\,{\nu}\,t-2\,h )(( 3/2\,{S_{32}}\,{c}^{2}-1/2\,{
				S_{32}}) \sin( g ) +{C_{32}} c\,\cos ( g)  )),
		\end{array}
	\end{equation*}
	
	
	
	\begin{equation*}\label{26}
		\begin{array}{ll}
			\langle R_{C33,S33} \rangle  & = {\frac {45}{4}}\sqrt{\frac{7}{360}}\,\frac {e{n}^{2}{R}^{3}_{M}}{{ (1 -{e}^{2})^{5/2}a}} ( 1/2{S_{33}}\,
			( c+1)( {c}^{2}-2/3\,c+1 ) \sin( {\nu}\,t-g-h) \\[1.5ex]
			\qquad & + ( c-1)( {C_{33}}\,\cos( g)( c+1) \cos( 3\,{\nu}\,t-3\,h
			) +{C_{33}}\,\sin ( g) c ( c+1) \sin ( 3\,{\nu}\,t-3\,h) \\[1.5ex]
			\qquad & - 1/2\, ( {c}^{2}+2/3\,c+1) {S_{33}}\,\sin ( {\nu}\,t+g-h))
			).
		\end{array}
	\end{equation*}
	
	The relation $\mu= n^{2}a^{3} = Gm$ is used, where $m$ is the mass of the central body and $G$ is the universal gravitational constant. In addition, $\nu t$ represents the angular velocity of the Moon orbit at time $t$. Finally, $a$ is the semimajor axis, $e$ represents the eccentricity, $i$ gives the inclination, $g$ is the argument of the perilune, $h$ is the longitude of the ascending node and, $n$ represents the mean motion.
	
	
	
	%%=============================================%%
	%% For submissions to Nature Portfolio Journals %%
	%% please use the heading ``Extended Data''.   %%
	%%=============================================%%
	
	%%=============================================================%%
	%% Sample for another appendix section			       %%
	%%=============================================================%%
	
	%% \section{Example of another appendix section}\label{secA2}%
	%% Appendices may be used for helpful, supporting or essential material that would otherwise 
	%% clutter, break up or be distracting to the text. Appendices can consist of sections, figures, 
	%% tables and equations etc.
	
\end{appendices}

%%===========================================================================================%%
%% If you are submitting to one of the Nature Portfolio journals, using the eJP submission   %%
%% system, please include the references within the manuscript file itself. You may do this  %%
%% by copying the reference list from your .bbl file, paste it into the main manuscript .tex %%
%% file, and delete the associated \verb+\bibliography+ commands.                            %%
%%===========================================================================================%%

%\bibliography{sn-bibliography}% common bib file
%% if required, the content of .bbl file can be included here once bbl is generated
%%\input sn-article.bbl


\begin{thebibliography}{}
	
	\bibitem[Abad et al.(2009)]{Abad} Abad, A., Elipe, A., \& Tresaco, E.: Analytical Model to Find Frozen Orbits for a Lunar Orbiter. J. Guid. Control Dyn., \textbf{32}(3), 888-898 (2009). https://doi.org/10.2514/1.38350
	
	\bibitem[Alessi et al.(2018)]{Alessi} Alessi, E. M., Schettino, G., Rossi, A., \& Valsecchi, G. B.: Solar radiation pressure resonances in Low Earth Orbits. MNRAS, \textbf{473}, 2407-2414 (2018). https://doi.org/10.1093/mnras/stx2507
    
    \bibitem[Carvalho et al.(2009)]{Carvalho6} Carvalho, J. P. S., Vilhena de Moraes, R., \& Prado, A. F. B. A.: Nonsphericity of the Moon and Near Sun-Synchronous Polar Lunar Orbits. Math. Probl. Eng., \textbf{2009}(1), Article ID 740460 (2009). https://doi.org/10.1155/2009/740460.

    \bibitem[Carvalho et al.(2010)]{Carvalho1} Carvalho, J. P. S., Vilhena de Moraes, R., \& Prado, A. F. B. A.: Some orbital characteristics of lunar artificial satellites. Celest. Mech. Dyn. Astron., \textbf{108}, 371-388 (2010). https://doi.org/10.1007/s10569-010-9310-6

    \bibitem[Carvalho et al.(2022a)]{Carvalho} Carvalho, J. P. S., Vilhena de Moraes, R., \& Prado, A. F. B. A.: Analysis of the orbital evolution of space debris using a solar sail and natural forces. Adv. Space Res., \textbf{70}, 125-143 (2022). https://doi.org/10.1016/j.asr.2022.04.014

    \bibitem[Carvalho et al.(2022b)]{Carvalho3} Carvalho, J. P. S., Yokoyama, T., \& Mour\~ao, D. C.: Single-averaged model for analysis of frozen orbits around planets and moons. Celest. Mech. Dyn. Astron., \textbf{134}(4), 35 (2022). https://doi.org/10.1007/s10569-022-10092-6

    \bibitem[Carvalho et al.(2023)]{Carvalho7} Carvalho, J. P. S., Tresaco, E., Elipe, A., \& Vilhena de Moraes, R.: Space Debris: A Mitigation Approach Based on Natural Perturbations and Resonances (2023). In: International Orbital Debris Conference. Available at: \url{https://www.hou.usra.edu/meetings/orbitaldebris2023/pdf/6084.pdf}.

    \bibitem[Coffey et al. (1994)]{Coffey} Coffey, S. L., Deprit, A. \& Deprit, E.: Frozen orbits for satellites close to an Earth-like planet. Celest. Mech Dyn Astron., \textbf{59}, 37–72 (1994). https://doi.org/10.1007/BF00691970

    \bibitem[De Melo et al.(2005)]{Melo} De Melo, C. F., Winter, O. C., \& Vieira Neto, E.: Numerical study of low-cost alternative orbits around the Moon. Adv. Space Res., \textbf{36}(3), 552-560 (2005). https://doi.org/10.1016/j.asr.2005.07.052.

    \bibitem[De Saedeleer(2006)]{Saedeleer} De Saedeleer B.: Analytical theory of lunar artificial satellite with third body perturbations. Celest. Mech. Dyn. Astron., \textbf{95}, 407-423 (2006). https://doi.org/10.1007/s10569-006-9029-6

    \bibitem[De Saedeleer \& Henrard (2006)]{Henrard} De Saedeleer, B., \& Henrard, J.: The combined effect of $J_2$ and $C_{22}$ on the critical inclination of a lunar orbiter. Adv. Space Res., \textbf{37}, 80-87 (2006). https://doi.org/10.1016/j.asr.2005.06.052

    \bibitem[Gomes \& Domingos (2016)]{Gomes} Gomes, V. M., \& De C\'assia Domingos, R.: Studying the lifetime of orbits around Moons in elliptic motion. Comp. Appl. Math., \textbf{35}, 653-661 (2016). https://doi.org/10.1007/s40314-015-0258-8

    \bibitem[Elipe \& Lara(2003)]{Elipe} Elipe, A., \& Lara, M.: Frozen Orbits About the Moon. J. Guid. Control Dyn., \textbf{26}(2), 238-243 (2003). https://doi.org/10.2514/2.5064

    \bibitem[El-Salam \& El-Bar(2016)]{ElSalam} El-Salam, F.A., \& El-Bar, S.E.: Families of frozen orbits of lunar artificial satellites. Appl. Math. Model., \textbf{40}(23-24), 9739-9753 (2016). https://doi.org/10.1016/j.apm.2016.06.036

    \bibitem[Folta \& Quinn(2006)]{Folta} Folta, D., \& Quinn, D.: Lunar Frozen Orbits. In: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, (2006). https://doi.org/10.2514/6.2006-6749

    \bibitem[Giacaglia et al.(1970)]{Giacaglia} Giacaglia, G.E.O., Murphy, J.P. \& Felsentreger, T.L.: A semi-analytic theory for the motion of a lunar satellite. Celest. Mech., \textbf{3}, 3–66 (1970). https://doi.org/10.1007/BF01230432

    \bibitem[Gon\c{c}alves et al.(2016)]{Liana} Gon\c{c}alves, L.D., Rocco, E.M., Vilhena de Moraes, R., \& Prado, A.F.B.A.: Effects of the individual terms of the lunar potential in the motion of satellites around the Moon. Int. J. Theor. Appl. Mech., \textbf{1}, 20-29 (2016). 

    \bibitem[Gupta \& Sharma(2011)]{Sharma} Gupta, S., \& Sharma, R.K.: Effect of Altitude, Right Ascension of Ascending Node and Inclination on Lifetime of Circular Lunar Orbits. Int. J. Astron. Astrophys, \textbf{1}, 155-163 (2011). https://doi.org/10.4236/ijaa.2011.13020   

    \bibitem[Kikuchi et al. (2021)]{Kikuchi} Kikuchi, S., Oki, Y., \& Tsuda, Y.: Frozen Orbits Under Radiation Pressure and Zonal Gravity Perturbations. Journal of Guidance, Control, and Dynamics, \textbf{44}(11), 1924-1946 (2021). https://doi.org/10.2514/1.G005564

    \bibitem[Konopliv et al.(1993)]{Konopliv} Konopliv, A.S., Sjogren, W., Wimberly, R.N., Cook, R.A. \& Vijayaraghavan, A.: A high resolution lunar gravity field and predicted orbit behavior (1993). In: AAS/AIAA Astrodynamics Specialist Conference, Paper AAS 93-622.

    \bibitem[Khattab et al.(2020)]{Khattab} Khattab, E.H., Radwan, M., \& Rahoma, W.A.: Frozen Orbits Construction for a Lunar Solar Sail. J. Astron. Space Sci., \textbf{37}(1), 1-9 (2020). https://doi.org/10.5140/JASS.2020.37.1.1

    \bibitem[Kubo-oka(1999)]{Kubo-oka} Kubo-oka, T.: Long-term effect of solar radiation on the orbit of an octagonal satellite orbiting around the Moon. Adv. Space Res., \textbf{23}(4), 727-731 (1999). https://doi.org/10.1016/S0273-1177(99)00146-5
 
	\bibitem[Kubo-oka \& Sengoku(1999)]{Kubo} Kubo-oka, T., \& Sengoku, A.: A. Solar radiation pressure model for the relay satellite of SELENE. Earth Planet Sp., \textbf{51}, 979–986 (1999). https://doi.org/10.1186/BF03351568

    \bibitem[Kuga et al.(2011)]{Kuga} Kuga, H.K., Carrara,V., \& Rao, K.R.: Sat\'elites artificiais - movimento orbital (2011). Instituto Nacional de Pesquisas Espaciais- INPE.

    
    \bibitem[Liu et al.(2022)]{Liu} Liu, J., Xu, B., Li, C., \& Li, M.: Lifetime Extension of Ultra Low-Altitude Lunar Spacecraft with Low-Thrust Propulsion System. Aerospace, \textbf{9}(6), 305 (2022). https://doi.org/10.3390/aerospace9060305

    \bibitem[Sirwah et al. (2020)]{Sirwah} Sirwah, M.A., Tarek, D., Radwan, M., \& Ibrahim, A.H.: A study of the moderate altitude frozen orbits around the Moon. Results Phys., \textbf{17}, 103 - 148 (2020). https://doi.org/10.1016/j.rinp.2020.103148 

    \bibitem[McInnes(1999)]{McInnes} McInnes, C.R.: Solar Sailing: Technology, Dynamics and Mission Applications. Springer-Verlag. (1999).

    \bibitem[Montenbruck et al.(2000)]{Montenbruck} Montenbruck, O., Gill, E., \& Mahooti, M.: Satellite Orbits Models Methods Applications, 1st edn. Springer-Verlag. (2000).

    \bibitem[Nie \& Gurfil(2018)]{Nie} Nie, T., \& Gurfil, P.: Lunar frozen orbits revisited. Celest Mech Dyn Astr., \textbf{130}, 61 (2018). https://doi.org/10.1007/s10569-018-9858-0

	
    
	\bibitem[Nie et al.(2025)]{Ren} Nie, T., Li, J., Zhang, S., Ren, J., \& Xu, R.: Analytical propagation for third-body perturbed orbits. Chin. J. Aeronaut., \textbf{39}, 4 (2025). https://doi.org/10.1016/j.cja.2025.103865

    \bibitem[Nie et al.(2026)]{Nie2026} Nie, T., Que, Z., Zhang, S., Ren, J., Xu, R.: Long-term frozen repeat orbits with large eccentricity under complex perturbations, Adv. Space Res., \textbf{77}(4), 4909-4927 (2026). https://doi.org/10.1016/j.asr.2025.11.106.
    
    \bibitem[Nie et al.(2024)]{Nie2024} Nie, T., Zhang, S., Yi, T., \& Ren, J.: Analytical Solution for Third-Body Perturbations with Double Averaging. J. Guid. Control Dyn., \textbf{48}(2), 384-392 (2024). https://doi.org/10.2514/1.G008338

    \bibitem[Oesterwinter(1970)]{Oesterwinter} Oesterwinter, C.: The Motion of a Lunar Satellite. Celest. Mech., \textbf{1}(3-4), 368-436 (1970). https://doi.org/10.1007/BF01231142.

    \bibitem[Pontani et al.(2018)]{Pontani} Pontani, M., Di Roberto, R., \& Graziani, F.: Lunar orbit dynamics and maneuvers for Lunisat missions. Acta Astronaut., \textbf{149}, 111-122 (2018).  https://doi.org/10.1016/j.actaastro.2018.05.015

    \bibitem[Popandopulo et al.(2022)]{Popandopulo} Popandopulo, N.A., Aleksandrova, A.G., Tomilova, I.V., Avdyushev, V.A., \& Bordovitsyna, T.V.: Numerical modelling of the dynamics of artificial satellites of the Moon. Sol. Syst. Res., \textbf{56}(4), 252–270 (2022). https://doi.org/10.31857/S0320930X22040077.

    \bibitem[Roncoli(2005)]{Roncoli} Roncoli, R.B.: Lunar Constants and Models Document. JPL Technical Document D-32296 (2005).

    \bibitem[Sirwah et al.(2020)]{Magdy} Sirwah, M.A., Tarek, D., Radwan, M., \& Ibrahim, A.H.: A study of the moderate altitude frozen orbits around the Moon. Results Phys., \textbf{17}, 103-148 (2020) https://doi.org/10.1016/j.rinp.2020.103148.

    \bibitem[Tresaco et al.(2016)]{Tresacob} Tresaco, E., Elipe, A., \& Carvalho, J.P.S.:  Frozen orbits for a solar sail around Mercury. J. Guid. Control Dyn., \textbf{39}(7), 1659-1666 (2016). https://doi.org/10.2514/1.G001510.

    \bibitem[Tresaco et al.(2018)]{Tresaco} Tresaco, E., Carvalho, J.P.S., Elipe, A., Prado, A.F.B.A., \& Vilhena de Moraes, R.: Averaged model to study long‑term dynamics of a probe about Mercury. Celest. Mech. Dyn. Astron., \textbf{130}, 9 (2018). https://doi.org/10.1007/s10569-017-9801-9

    \bibitem[Tzirti et al.(2014)]{Tzirti} Tzirti, S., Noullez, A., \& Tsiganis, K.: Secular dynamics of a lunar orbiter: a global exploration using Prony’s frequency analysis. Celest. Mech. Dyn. Astr., \textbf{118}, 379–397 (2014). https://doi.org/10.1007/s10569-014-9540-0
    
    \bibitem[Vinti(1971)]{Vinti} Vinti, J. P. (1971).: Representation of the Earth’s Gravitational Potential. Celest. Mech., \textbf{4}(3-4), 348 - 367 (1971). https://doi.org/10.1007/BF01231397

    \bibitem[Xu(2025)]{Xu2025} Xu, X.: Analysis of Near-Polar and Near-Circular Periodic Orbits Around the Moon with $J_2$, $C_{22}$ and Third-Body Perturbations. Symmetry 2025, \textbf{17}, 630 (2025). https://doi.org/10.3390/sym17050630
    
\end{thebibliography}


\end{document}

