Abstract Artificial satellite theories for the Moon differ in important aspects from the similar theories for Earth because of more uneven nature of lunar gravity, strong third-body perturbations, and no atmospheric drag. The existing lunar satellite theories include only a few gravity spherical harmonics along with Hill’s approximation of the third-body disturbing function. Despite such simplifications, they provide semi-analytical solutions and require numerical propagation, thus limiting their applications. In this work, a fully analytic theory for lunar satellites is developed using generalized formulae that enable inclusion of spherical harmonics up to an arbitrary degree and order. The complete gravity and third-body disturbing functions are together considered as the dominant perturbation to the Keplerian Hamiltonian, which is then fully normalized up to second order by constructing three canonical transformations for averaging out the short-, medium-, and long-period terms. The gravity harmonics are fully treated in Part-I of this paper with the short- and medium-period terms developed up to first order and the long-period and secular terms up to second order. These second-order terms capture the coupling effects of the gravity harmonics. For the third-body perturbation, the secular terms are developed up to first order here. The periodic effects (including those originating from the coupling of gravity and third-body terms) will be explicitly treated in Part-II of this paper. The significance of the coupling effects from the gravity harmonics is assessed by comparing the analytically-propagated ephemerides of low-altitude lunar satellites with the numerically-propagated results. 2 Bharat Mahajan 1 Introduction Lunar missions have gained renewed interest recently from many government and private space organizations. New lunar campaigns are being designed with evolution- ary architectures enabling long-term human presence on the Moon ( Condon et al. , 2020). Additional space infrastructure supporting communication and navigation on the lunar surface will inevitably be required for achieving this goal. The number of space assets in lunar orbits, therefore, are expected to increase in the near fu- ture and this will necessitate the availability of more accurate and efficient models for lunar orbit dynamics. While special perturbation methods are capable of accu- rate orbital state predictions, they offer limited insight into the long-term dynam- ics of artificial satellites. General perturbation methods, on the other hand, provide computationally-efficient techniques for long-term state and uncertainty propagation ( Hoots and Roehrich, 1980; Mahajan, 2018), orbit maintenance including guidance & control (Kneˇzevi´c and Milani , 1998; Alfriend et al. , 2009), lifetime estimation ( Leg- naro and Efthymiopoulos , 2024), and frozen orbit design ( Abad et al. , 2009; Lara, 2011; Nie and Gurfil , 2018) to name a few applications. The most significant perturbations affecting artificial satellites in lunar orbits are nonuniform gravity and third-body effects. These two perturbations have more pro- nounced effects in lunar orbits than Earth orbits. This is due to two reasons: more uneven mass distribution (read: mascons) of the Moon including less flattening of its shape and the lunar sphere of influence in the Earth-Moon system being much smaller than that of Earth. As a consequence of the former, multiple zonal, sectoral, and tesseral spherical harmonic coefficients in the lunar gravity models are of the same order of magnitude and they are only one order of magnitude smaller than that of the oblateness coefficient. For Earth, the latter is at least three order of magnitude larger than any other coefficient. Satellites in medium-to-high altitude lunar orbits also experience more significant third-body gravitational effects from Earth than vice versa. For constructing artificial satellite theories for the Moon, the most significant spherical harmonics must be together considered as the dominant perturbation to the Keplerian motion. However, this may cause difficulties by using the conventional approaches for separating the long-periodic effects from the equations of motion to enable analytic integration. Most theories for lunar satellites in the literature are semi- analytical as the long-periodic effects are not separated and they require numerical integration (see Kozai (1963), Giacaglia et al. (1970), De Saedeleer (2006), San- Juan et al. (2019), Efthymiopoulos et al. (2023) and references therein). Kneˇzevi´c and Milani (1995) obtained a first-order fully analytic theory without separating the long-periodic effects but only for the special Lunar satellite analytic theory with complete gravity and third-body perturbations 11 Fig. 1 Disturbing function magnitudes (scaled by the Keplerian Hamiltonian term) versus the lunar orbit altitude for the nonuniform gravity and third-body perturbations. The subscripts indicate degree and order of the spherical harmonics. Table 1 Mean values of the astronomical quantities used for computations. Constant V alue Description R/leftmoon1738 Radius (GRGM1200A) of the Moon µ /leftmoon4902.8001224453001 Gravitational parameter (GRGM1200A) of the Moon a⨁ 384400 Semimajor axis of the Earth’s orbit in the Moon-centered inertial frame µ ⨁ 398600.4415 Gravitational parameter of Earth e⨁ 0.0549 Eccentricity of the Earth’s orbit in the Moon-centered inertial frame a⊙ 149597871 Semimajor axis of the Sun’s orbit in the Moon-centered inertial frame µ ⊙ 1.32712440018 × 1011 Gravitational parameter of the Sun e⊙ 0.01671022 Eccentricity of the Sun’s orbit in the Moon-centered inertial frame where H0 = HKep, H1 = −R′ − R. The formal parameter ε is only used to separate the terms in accordance with their order and is set to unity in the final results. The canonical perturbation method used for normalizing this Hamiltonian up to second order is described in the next section.