Luna 3 (or Lunik 3 in Russian sources) was the first spacecraft to perform a flyby of the Moon. Launched in October 1959 on a translunar trajectory with large semi-major axis and eccentricity, it collided with the Earth in late March 1960. The short, 6-month dynamical lifetime has often been explained through an increase in eccentricity due to the Lidov-Kozai effect. However, the classical Lidov-Kozai solution is only valid in the limit of small semi-major axis ratio, a condition that is satisfied only for solar (but not for lunar) perturbations. We undertook a study of the dynamics of Luna 3 with the aim of assessing the principal mechanisms affecting its evolution. We analyze the Luna 3 trajectory by generating accurate osculating solutions, and by comparing them to integrations of singly- and doubly-averaged equations of motion in vectorial form. Lunar close encounters, which cannot be reproduced in an averaging approach, decisively affect the trajectory and break the doubly-averaged dynamics. Solar perturbations induce oscillations of intermediate period that affect the geometry of the close encounters and cause the singly-averaged and osculating inclinations to change quadrants (the orbital plane ``flips''). We find that the peculiar evolution of Luna 3 can only be explained by taking into account lunar close encounters and intermediate-period terms; such terms are averaged out in the Lidov-Kozai solution, which is not adequate to describe translunar or cislunar trajectories. Understanding the limits of the Lidov-Kozai solution is of particular significance for the motion of objects in the Earth-Moon environment and of exoplanetary systems.

We aim to provide satellite operators and researchers with an efficient means for evaluating and mitigating collision risk during the design process of mega-constellations. We first establish a baseline for evaluating various techniques for close-encounter prediction and collision-probability calculation (Hoots et al. 1984, Gronchi 2005, JeongAhn and Malhotra 2015) by carrying out brute-force numerical simulations and using a sequence of filters to greatly reduce the computational expense of the algorithm. Next, we estimate conjunction events in the orbital environment following the anticipated deployments of the OneWeb LEO and SpaceX Starlink mega-constellations. As a final step, we investigate Minimum Space Occupancy (MiSO) orbits (Bombardelli et al. 2018), a generalization of the well-known frozen orbits that account for the perturbed-Keplerian dynamics of the Earth-Moon-Sun-satellite system. We evaluate the ability of MiSO configurations of the proposed mega-constellations, as suggested by Bombardelli et al. 2018, to reduce the risk of endogenous (intra-constellation) collisions. The results indicate that the adoption of the MiSO orbital configuration can significantly reduce risk with nearly indistinguishable adjustments to the nominal orbital elements of the constellation satellites.

The description of the long-term dynamics of highly elliptic orbits under third-body perturbations may require an expansion of the disturbing function in series of the semi-major axes ratio up to higher orders. To avoid dealing with long series in trigonometric functions, we refer the motion to the apsidal frame and efficiently remove the short-period effects of this expansion in vectorial form up to an arbitrary order. We then provide the variation equations of the two fundamental vectors of the Keplerian motion by analogous vectorial recurrences, which are free from singularities and take a compact form useful for the numerical propagation of the flow in mean elements.

This paper is concerned with the comparison of semi-analytical and non-averaged propagation methods for Earth satellite orbits. We analyse the total integration error for semi-analytical methods and propose a novel decomposition into dynamical, model truncation, short-periodic, and numerical error components. The first three are attributable to distinct approximations required by the method of averaging, which fundamentally limit the attainable accuracy. In contrast, numerical error, the only component present in non-averaged methods, can be significantly mitigated by employing adaptive numerical algorithms and regularized formulations of the equations of motion. We present a collection of non-averaged methods based on the integration of existing regularized formulations of the equations of motion through an adaptive solver. We implemented the collection in the orbit propagation code THALASSA, which we make publicly available, and we compared the non-averaged methods to the semi-analytical method implemented in the orbit propagation tool STELA through numerical tests involving long-term propagations (on the order of decades) of LEO, GTO, and high-altitude HEO orbits. For the test cases considered, regularized non-averaged methods were found to be up to two times slower than semi-analytical for the LEO orbit, to have comparable speed for the GTO, and to be ten times as fast for the HEO (for the same accuracy). We show for the first time that efficient implementations of non-averaged regularized formulations of the equations of motion, and especially of non-singular element methods, are attractive candidates for the long-term study of high-altitude and highly elliptical Earth satellite orbits.

We have carried out a numerical investigation of the coupled gravitational and non-gravitational perturbations acting on Earth satellite orbits in an extensive grid, covering the whole circumterrestrial space, using an appropriately modified version of the SWIFT symplectic integrator, which is suitable for long-term (120 years) integrations of the non-averaged equations of motion. Hence, we characterize the long-term dynamics and the phase-space structure of the Earth-orbiter environment, starting from low altitudes (400 km) and going up to the GEO region and beyond. This investigation was done in the framework of the EC-funded "ReDSHIFT" project, with the purpose of enabling the definition of passive debris removal strategies, based on the use of physical mechanisms inherent in the complex dynamics of the problem (i.e., resonances). Accordingly, the complicated interactions among resonances, generated by different perturbing forces (i.e., lunisolar gravity, solar radiation pressure, tesseral harmonics in the geopotential) are accurately depicted in our results, where we can identify the regions of phase space where the motion is regular and long-term stable and regions for which eccentricity growth and even instability due to chaotic behavior can emerge. The results are presented in an "atlas" of dynamical stability maps for different orbital zones, with a particular focus on the (drag-free) range of semimajor axes, where the perturbing effects of the Earth's oblateness and lunisolar gravity are of comparable order. In some regions, the overlapping of the predominant lunisolar secular and semi-secular resonances furnish a number of interesting disposal hatches at moderate to low eccentricity orbits. All computations were repeated for an increased area-to-mass ratio, simulating the case of a satellite equipped with an on-board, area-augmenting device.

In this paper, we study the long-term dynamical evolution of highly-elliptical orbits (HEOs) in the medium-Earth orbit (MEO) region around the Earth. The real population consists primarily of Geosynchronous Transfer Orbits (GTOs), launched at specific inclinations, Molniya-type satellites and related debris. We performed a suite of long-term numerical integrations (up to 200 years) within a realistic dynamical model, aimed primarily at recording the dynamical lifetime of such orbits (defined as the time needed for atmospheric reentry) and understanding its dependence on initial conditions and other parameters, such as the area-to-mass ratio (A/m). Our results are presented in the form of 2-D lifetime maps, for different values of inclination, A/m, and drag coefficient. We find that the majority of small debris (> 70%, depending on the inclination) can naturally reenter within 25-90 years, but these numbers are significantly less optimistic for large debris (e.g., upper stages), with the notable exception of those launched from high latitude (Baikonur). We estimate the reentry probability and mean dynamical lifetime for different classes of GTOs and we find that both quantities depend primarily and strongly on initial perigee altitude. Atmospheric drag and higher A/m values extend the reentry zones, especially at low inclinations. For high inclinations, this dependence is weakened, as the primary mechanisms leading to reentry are overlapping lunisolar resonances. This study forms part of the EC-funded (H2020) "ReDSHIFT" project.

The Medium Earth Orbit (MEO) region hosts satellites for navigation, communication, and geodetic/space environmental science, among which are the Global Navigation Satellites Systems (GNSS). Safe and efficient removal of debris from MEO is problematic due to the high cost for maneuvers needed to directly reach the Earth (reentry orbits) and the relatively crowded GNSS neighborhood (graveyard orbits). Recent studies have highlighted the complicated secular dynamics in the MEO region, but also the possibility of exploiting these dynamics, for designing removal strategies. In this paper, we present our numerical exploration of the long-term dynamics in MEO, performed with the purpose of unveiling the set of reentry and graveyard solutions that could be reached with maneuvers of reasonable DV cost. We simulated the dynamics over 120-200 years for an extended grid of millions of fictitious MEO satellites that covered all inclinations from 0 to 90deg, using non-averaged equations of motion and a suitable dynamical model that accounted for the principal geopotential terms, 3rd-body perturbations and solar radiation pressure (SRP). We found a sizeable set of usable solutions with reentry times that exceed ~40years, mainly around three specific inclination values: 46deg, 56deg, and 68deg; a result compatible with our understanding of MEO secular dynamics. For DV <= 300 m/s (i.e., achieved if you start from a typical GNSS orbit and target a disposal orbit with e<0.3), reentry times from GNSS altitudes exceed ~70 years, while low-cost (DV ~= 5-35 m/s) graveyard orbits, stable for at lest 200 years, are found for eccentricities up to e~0.018. This investigation was carried out in the framework of the EC-funded "ReDSHIFT" project.

The slow deformation of terrestrial orbits in the medium range, subject to lunisolar resonances, is well approximated by a family of Hamiltonian flow with $2.5$ degree-of-freedom. The action variables of the system may experience chaotic variations and large drift that we may quantify. Using variational chaos indicators, we compute high-resolution portraits of the action space. Such refined meshes allow to reveal the existence of tori and structures filling chaotic regions. Our elaborate computations allow us to isolate precise initial conditions near specific zones of interest and study their asymptotic behaviour in time. Borrowing classical techniques of phase- space visualisation, we highlight how the drift is mediated by the complement of the numerically detected KAM tori.

We consider Earth satellite orbits in the range of semi-major axes where the perturbing effects of Earth's oblateness and lunisolar gravity are of comparable order. This range covers the medium-Earth orbits (MEO) of the Global Navigation Satellite Systems and the geosynchronous orbits (GEO) of the communication satellites. We recall a secular and quadrupolar model, based on the Milankovitch vector formulation of perturbation theory, which governs the long-term orbital evolution subject to the predominant gravitational interactions. We study the global dynamics of this two-and-a-half degrees-of-freedom Hamiltonian system by means of the fast Lyapunov indicator (FLI), used in a statistical sense. Specifically, we characterize the degree of chaoticity of the action space using angle-averaged normalized FLI maps, thereby overcoming the angle dependencies of the conventional stability maps. Emphasis is placed upon the phase-space structures near secular resonances, which are of first importance to the space debris community. We confirm and quantify the transition from order to chaos in MEO, stemming from the critical inclinations, and find that highly inclined GEO orbits are particularly unstable. Despite their reputed normality, Earth satellite orbits can possess an extraordinarily rich spectrum of dynamical behaviors, and, from a mathematical perspective, have all the complications that make them very interesting candidates for testing the modern tools of chaos theory.

We provide a detailed derivation of the analytical expansion of the lunar and solar disturbing functions. Although there exist several papers on this topic, many derivations contain mistakes in the final expansion or rather (just) in the proof, thereby necessitating a recasting and correction of the original derivation. In this work, we provide a self-consistent and definite form of the lunisolar expansion. We start with Kaula's expansion of the disturbing function in terms of the equatorial elements of both the perturbed and perturbing bodies. Then we give a detailed proof of Lane's expansion, in which the elements of the Moon are referred to the ecliptic plane. Using this approach the inclination of the Moon becomes nearly constant, while the argument of perihelion, the longitude of the ascending node, and the mean anomaly vary linearly with time. We make a comparison between the different expansions and we profit from such discussion to point out some mistakes in the existing literature, which might compromise the correctness of the results. As an application, we analyze the long--term motion of the highly elliptical and critically--inclined Molniya orbits subject to quadrupolar gravitational interactions. The analytical expansions presented herein are very powerful with respect to dynamical studies based on Cartesian equations, because they quickly allow for a more holistic and intuitively understandable picture of the dynamics.

The navigation satellite constellations in medium-Earth orbit exist in a background of third-body secular resonances stemming from the perturbing gravitational effects of the Moon and the Sun. The resulting chaotic motions, emanating from the overlapping of neighboring resonant harmonics, induce especially strong perturbations on the orbital eccentricity, which can be transported to large values, thereby increasing the collision risk to the constellations and possibly leading to a proliferation of space debris. We show here that this transport is of a diffusive nature and we present representative diffusion maps that are useful in obtaining a global comprehension of the dynamical structure of the navigation satellite orbits.

Recent studies have shown that the medium-Earth orbit (MEO) region of the Global Navigation Satellite Systems is permeated by a devious network of lunisolar secular resonances, which can interact to produce chaotic and diffusive motions. The precarious state of the four navigation constellations, perched on the threshold of instability, makes it understandable why all past efforts to define stable graveyard orbits, especially in the case of Galileo, were bound to fail; the region is far too complex to allow of an adoption of the simple geosynchronous disposal strategy. We retrace one such recent attempt, funded by ESA's General Studies Programme in the frame of the GreenOPS initiative, that uses a systematic parametric approach and the straightforward maximum-eccentricity method to identify long-term stable regions, suitable for graveyards, as well as large-scale excursions in eccentricity, which can be used for post-mission deorbiting of constellation satellites. We then apply our new results on the stunningly rich dynamical structure of the MEO region toward the analysis of these disposal strategies for Galileo, and discuss the practical implications of resonances and chaos in this regime. We outline how the identification of the hyperbolic and elliptic fixed points of the resonances near Galileo can lead to explicit criteria for defining optimal disposal strategies.

It has long been suspected that the Global Navigation Satellite Systems exist in a background of complex resonances and chaotic motion; yet, the precise dynamical character of these phenomena remains elusive. Recent studies have shown that the occurrence and nature of the resonances driving these dynamics depend chiefly on the frequencies of nodal and apsidal precession and the rate of regression of the Moon's nodes. Woven throughout the inclination and eccentricity phase space is an exceedingly complicated web-like structure of lunisolar secular resonances, which become particularly dense near the inclinations of the navigation satellite orbits. A clear picture of the physical significance of these resonances is of considerable practical interest for the design of disposal strategies for the four constellations. Here we present analytical and semi-analytical models that accurately reflect the true nature of the resonant interactions, and trace the topological organization of the manifolds on which the chaotic motions take place. We present an atlas of FLI stability maps, showing the extent of the chaotic regions of the phase space, computed through a hierarchy of more realistic, and more complicated, models, and compare the chaotic zones in these charts with the analytical estimation of the width of the chaotic layers from the heuristic Chirikov resonance-overlap criterion. As the semi-major axis of the satellite is receding, we observe a transition from stable Nekhoroshev-like structures at three Earth radii, where regular orbits dominate, to a Chirikov regime where resonances overlap at five Earth radii. From a numerical estimation of the Lyapunov times, we find that many of the inclined, nearly circular orbits of the navigation satellites are strongly chaotic and that their dynamics are unpredictable on decadal timescales.

Numerical simulations carried out over the past decade suggest that the orbits of the Global Navigation Satellite Systems are unstable, resulting in an apparent chaotic growth of the eccentricity. Here we show that the irregular and haphazard character of these orbits reflects a similar irregularity in the orbits of many celestial bodies in our Solar System. We find that secular resonances, involving linear combinations of the frequencies of nodal and apsidal precession and the rate of regression of lunar nodes, occur in profusion so that the phase space is threaded by a devious stochastic web. As in all cases in the Solar System, chaos ensues where resonances overlap. These results may be significant for the analysis of disposal strategies for the four constellations in this precarious region of space.