Frozen or quasi-frozen initial conditions are indispensable in order to provide stable lunar orbits. This is a mainly a consequence of the complex gravitational field of the Moon combined to its slow rotation compared to the Earth one. Additionally, while the Earth gravitational environment is mainly dominated by a few terms, Lunar orbits are influenced by a high degree and order of zonal and tesseral harmonics and third body effects, which makes them very sensitive to initial conditions and greatly limits their stability. Frozen orbits around the Earth can be easily obtained using, for example, the theory developed by Graham Cook in the 60s or straight from the 1959 Kozai's solution of the main problem (with the first harmonics) together with Lyddane regularization. On the other hand, when considering the lunar environment the design of frozen orbits gets considerably more complicated, greatly challenging the development of analytical methods. In this paper, an extensive analysis of low-altitude frozen orbits around the Moon is conducted using a combined approach including analytical tools and high-fidelity numerical propagation. More specifically, we employ an extended Kozai-Lyddane transformation by generalizing the work of Kozai 1959 to an arbitrary number of zonal harmonics (truncated to 50 for reasonable computational speed), which allows us to compute highly frozen solutions of a full-zonal model for generic inclination and semimajor axis. Because no restrictions are imposed on the orbit eccentricity, this approach not only generalizes Kozai's and Brouwer's work but also Kaula's work, which is only accurate for small eccentricities and can be considerably expensive from the numerical computation point of view when dealing with high order gravitational models.
After confirming the suitability of the developed analytical approach to obtain frozen orbit initial conditions in a full-zonal model we move to a high-fidelity model including tesseral harmonics up to 150x150, third body effects by the Sun and Earth and accounting for Lunar spin axis precession. A novel approach to evaluate the effective “temperature” of the resulting quasi-frozen orbit is adopted based on the recently proposed concept of “space occupancy” by one of the authors. Interestingly, low altitude lunar orbits resulting from full-zonal frozen conditions exhibits a lifetime of a few years in the high-fidelity model. On the other hand, higher eccentricity frozen orbits are less stable mainly owing to third-body effects, which suggests possible improvements of the proposed formulation in future works. The presented results highlight the effectiveness of the hybrid analytical and numerical approach for rapid analysis of highly stable quasi-frozen orbits across a wide parameter space.