Abstract
Satellite collisions or fragmentations generate a huge number of space debris; over time, the fragments might get dispersed, making it difficult to associate them to the configuration at breakup. In this work, we present a procedure to backtrace the debris, reconnecting them to their original configuration. To this end, we compute the proper elements, namely dynamical quantities which stay nearly constant over time. While the osculating elements might spread and lose connection with the values at breakup, the proper elements, which have been already successfully used to identify asteroid families, retain the dynamical features of the original configuration. We show the efficacy of the procedure, based on a hierarchical implementation of perturbation theory, by analyzing the following four different case studies associated to satellites that underwent a catastrophic event: Ariane 44lp, Atlas V Centaur, CZ3, Titan IIIc Transtage. The link between (initial and final) osculating and proper elements is evaluated through tools of statistical data analysis. The results show that proper elements allow one to reconnect the fragments to their parent body.
Introduction
Since the launch of Sputnik 1 in 1957, thousands of satellites have been deployed in orbit around the Earth. Explosions or collisions of satellites generated millions of space debris of various sizes^{1}, currently traveling at different altitudes: rocket stages, fragments from disintegrations, bolts, paint flakes, electronic parts, etc. Chain reactions triggered by catastrophic events involving satellites might increase the hazard of (human and robotic) space activities. A single breakup event generates a cloud of debris which scatters around, sometimes reaching great distances after a relatively short time. Once the fragments are dispersed, it is often difficult to trace them back; hence, a question of paramount importance is to connect the debris to their parent satellite. In this work we propose a method that allows us to link the fragments, after a certain interval of time, to the configuration of the debris soon after the initial catastrophic event. This result contributes to address a timely problem since, in case of a collision between two satellites or an explosion of a single satellite, it is certainly important to know the parent bodies that generated the space debris. The implications are wide and range from space sustainability to space law.
To study a specific breakup event, we introduce a suitable model (based on the Hamiltonian formalism) to describe the dynamics of each fragment. The model is composed of the sum of the Keplerian attraction, the effect of the geopotential, the gravitational influence of Sun and Moon. Then, we implement perturbation theory to construct a sequence of canonical transformations providing, for each debris, approximate integrals of motion called proper elements, namely quantities that stay nearly constant over time. Each fragment is characterized by a set of six orbital elements, namely semimajor axis a, eccentricity e, inclination i, mean anomaly M, argument of perigee \(\omega\) and longitude of the ascending node \(\Omega\). Starting from their initial values, we compute the orbital elements of each fragment after a given interval of time, to which we refer as the final osculating elements. Then, we compute the proper elements associated to the final osculating elements, and we compare them either with the initial elements and with the corresponding proper elements at the initial time. The comparison gives the desired information: while the final osculating elements might spread far away from the initial values, the (initial/final) proper elements stay almost constant and retain the original features of the cloud of fragments^{2}. A striking use of proper elements was already proposed to group asteroids, inspired by the pioneer work^{3} of Hirayama in 1918 and continued by many other authors^{4,5}. The analytical computation of proper elements allowed to group asteroids in families, possibly leading to the conjecture that such asteroids might be fragments of an ancestor parent body. Knežević and Milani^{6} introduced also the synthetic proper elements based upon a numerical integration, a digital filtering of the shortperiod terms and a Fourier analysis.
Motivated by the successful results on asteroids, we propose to group and reconnect space debris through the computation of the proper elements associated to fragments generated by a satellite breakup event^{7,8,9}. The procedure we are going to describe, requires the introduction of a realistic model which depends on quantities varying on different time scales; hence we need a suitable hierarchical set of transformations of coordinates, called normal forms, aimed at constructing the proper elements, whose relation with the initial elements is analyzed through statistical methods. We will consider four sample cases associated to the breakup events of the satellites Ariane 44lp, Atlas V Centaur, CZ3, Titan IIIc Transtage. Using statistical data analysis, we show the effectiveness of the use of proper elements in reconnecting the fragments to their parent body. To reconnect the debris to a parent body, we backpropagate the debris for a given time and compare the osculating or proper elements at the initial time and at the backpropagated time. The effectiveness of the method has been shown in the specific example of Titan III Transtage. We finally provide an example in which one can distinguish between proper elements associated to nearby breakup events.
This work is organized as follows. After introducing the model, we describe the procedure to compute the proper elements through normal form theory. Then we investigate the test cases by computing osculating and proper elements, and by analyzing the results through histograms, KolmogorovSmirnov test, Variance Equivalence test and Pearson correlation coefficient. We end up with some conclusions and perspectives.
The model
For the present work, our case studies will be located at altitudes between 15000 and 25000 km, all of them well above the Earth’s atmosphere. At those altitudes a celestial object is subject to different forces that we describe through a Hamiltonian function composed by the following parts: the attraction of the Earth (that we split as the sum of the Keplerian part \(H_{Kep}\) and the potential \(H_E\) generated by the Earth’s nonspherical shape), and the gravitational influence of the Moon \(H_M\) and the Sun \(H_S\) (both assumed to be point masses). The overall Hamiltonian
depends upon the orbital elements of the debris, Moon and Sun, and on the sidereal time describing the rotation of the Earth^{10}.
We are aware that a realistic model should include also the effect of the solar radiation pressure^{11} (SRP). However, we decided not to consider SRP for two main reasons: (i) the work^{2} provides some experiments on synthetic space debris (namely obtained through a simulator of breakup events), using a model that includes SRP; however, the results show that at intermediate altitudes the computation of the proper elements is not much affected by SRP, at least for objects with an areatomass ratio lower than 0.74; (ii) there does not exist a public catalogue that provides information about the areatomass ratio of real space debris, thus preventing reliable experiments on real cases.
The Keplerian and geopotential Hamiltonians
Expressing the Hamiltonian in terms of the orbital elements, the Keplerian part \(H_{Kep}\) is given by
where G is the gravitational constant and \(M_E\) is the mass of the Earth.
The contribution \(H_E\) due to the Earth’s nonspherical shape is computed as follows^{12,13}: we expand the geopotential in spherical harmonics, then we average over the fast variables (namely the mean anomaly of the debris and the sidereal time), and finally we limit the expansion of the secular part of the geopotential to the greatest spherical harmonic coefficients, usually denoted as \(J_2\) and \(J_3\). The resulting Hamiltonian takes the form:
where \(R_E\) is the Earth’s radius (equal to 6378.1363 km).
Moon and Sun Hamiltonians
The Hamiltonians of the Moon \(H_M\) and of the Sun \(H_S\) are expanded in powers of the ratios of the semimajor axes of the debris, and of the Moon and Sun, respectively \(a_M\) and \(a_S\), as well as in powers of the eccentricity and of the cosine of the inclination. The resulting expansions are truncated to a low order (typically the second one). By \((a_b,e_b,i_b,M_b,\omega _b,\Omega _b)\) we denote the orbital elements of the third body perturber, where \(b=M\) and \(b=S\) correspond to the Moon and Sun, respectively.
The expansion of the Moon’s Hamiltonian in terms of the orbital elements of the Moon and the debris is given below; we underline that in applications we will consider the expansion of \(H_M\) in (1) to \(l=2\):
where \(m_M\) is the mass of the Moon, \(y_s = 0\), if (s mod 2)=0, \(y_s = \frac{1}{2}\), if (s mod 2)= 1, \(t = (l1)\) mod 2, and
The functions \(F_{lmp}(i)\), \(F_{lsq}(i_M)\) and \(G_{lqr}(e_M)\) can be found, e.g., in^{10,12}; \(H_{lpj}(e)\) are the Hansen coefficients, while the terms \(U_l^{m,s}\) are given by
where \(\epsilon = 23^o26'21.406''\) is the Earth’s obliquity.
The Hamiltonian due to the Sun depends on the orbital elements of the Sun and the debris. The expansion of \(H_S\) is given below and, again, we will consider the expansion to \(l=2\):
where \(m_S\) is the mass of the Sun and
Validation of the model
Although representing an approximation of the full model, the Hamiltonian H in (1) provides an accurate description of the dynamics, as shown by the comparison in Fig. 1 between the integration of Hamilton’s equations and the Cartesian equations of motion^{13}.
Besides depending on the orbital elements of the debris, the Hamiltonian depends also upon the orbital elements of Moon and Sun. For our purposes, it is essential to stress that the debris, Moon and Sun move on different timescales, since the angular variables describing their respective motions vary with rates of the order of days (for the debris), months (for the Moon), years (for the Sun), see Table 1. As a consequence, the respective angular variables of debris, Moon and Sun can be ordered hierarchically as fast, semifast and slow. The fast angles are indeed the mean anomaly of the debris and the sidereal time accounting for the rotation of the Earth; we report in Fig. 1 also the integration obtained using the Hamiltonian, doubly averaged with respect to such fast angles.
Normal form and proper elements
We briefly recall the basics of normal form theory^{14}, which is at the basis of the computation of the proper elements. We consider a Hamiltonian of the form
where \(({{\underline{I}}},{{{\underline{\varphi }}}})\) are actionangle variables with \(({{\underline{I}}},{{{\underline{\varphi }}}})\in B\times {{{\mathbb {T}}}}^n\), where \(B\subset {{{\mathbb {R}}}}^n\) is an open set and n denotes the number of degrees of freedom. In (2) the function \({\mathscr {H}}_0({{\underline{I}}})\) is the integrable part, \(\varepsilon \in {{\mathbb {R}}}\) is a small parameter, \({\mathscr {H}}_1({{\underline{I}}},{{{\underline{\varphi }}}})\) is the perturbing function.
The normalization procedure consists in the definition of a suitable change of coordinates that transforms the Hamiltonian, so that it becomes integrable to orders of \(\varepsilon ^2\). The procedure can be iterated for some steps, but it is known that in general it is not converging^{15}.
We assume that the function \({\mathscr {H}}_1\) can be expanded in Fourier series as
where \(K\subseteq {{{\mathbb {Z}}}}^n\) and \(b_{{\underline{k}}}\) are functions with real coefficients. Let \(\chi\) be the generating function of the canonical transformation from the variables \(({{\underline{I}}},{{{\underline{\varphi }}}})\) to the new variables \(({\underline{I}}',{{{\underline{\varphi }}}}')\) given by
where the action of \(S_{\chi }^{\varepsilon }\) is defined by
with \(\{\cdot ,\cdot \}\) the Poisson bracket operator. We determine \(S_{\chi }^{\varepsilon }\) by requiring that the new Hamiltonian \({\mathscr {H}}^{(1)} = S_{\chi }^{\varepsilon }{\mathscr {H}}\) is transformed into
where \(Z_1={\mathscr {H}}_0+\varepsilon \overline{{\mathscr {H}}}_1\) is integrable (\(\overline{{\mathscr {H}}}_1\) is the average of \({\mathscr {H}}_1\)) and \({\mathscr {H}}_2\) is the remainder term of order \(\varepsilon ^2\). Inserting the change of coordinates in (2), one obtains the transformed Hamiltonian which takes the desired form (3) provided \(\chi\) satisfies the following normal form equation:
Expanding \(\chi\) in Fourier series, denoting the frequency by \({\underline{\omega }}_0= {{\partial {\mathscr {H}}_0}\over {\partial {{\underline{I}}}'}}\), one obtains that the generating function is given by the following formula, which is valid under the nonresonance assumption \({{\underline{k}}\cdot {{{\underline{\omega }}}}_0}\not =0\):
A higher order normal form is obtained by iterating the above procedure.
Recalling that the space debris model described above depends on fast, semifast and slow variables, we compute the normal form, taking advantage of the hierarchical structure of the coordinates associated to the debris, Moon and Sun. We first average the Hamiltonian over the fast (mean anomaly of the debris and sidereal time) and semifast (mean anomalies of Moon and Sun) angles. According to Hamilton’s equations, the rate of variation of the semimajor axis of the debris is given by the derivative of the Hamiltonian with respect to the mean anomaly; since we averaged over the mean anomaly, the semimajor axis is constant and becomes the first proper element, namely a quasiintegral of motion for the averaged approximated model. After averaging over the mean anomalies and the sidereal time, we endup with a Hamiltonian function with three degreesoffreedom in the extended phase space, since the Hamiltonian depends on time through the variation of the longitude of the ascending node of the Moon (see Table 1).
Next, we consider some reference values for the eccentricity and the inclination (namely the values of the fragments of the case study) and we expand the averaged Hamiltonian around such values. Then, we implement a canonical change of variables through a Lie series normalization, implemented through a Mathematica\(^{\copyright }\) program, that removes the dependence on the angles; this procedure provides two more proper elements associated with the eccentricity and the inclination. By making explicit all transformations^{2}, we end the procedure by backtransforming the change of variables to express the proper elements in the original coordinates.
In conclusion, the procedure leading to the computation of the proper elements can be summarized as follows^{2}.

1.
We consider the Hamiltonian including the contributions of the gravitational attractions of the Earth, Moon and Sun; we average with respect to the fast variables, in particular the mean anomaly M; hence, the semimajor axis is constant and becomes the first proper element.

2.
Since the longitude of the ascending node of the Moon \(\Omega _M\) depends on time, the Hamiltonian resulting from step 1 depends on \((e,i,\omega ,\Omega ,t)\); hence, we introduce the Hamiltonian in the extended phase space, so that it becomes autonomous, although depending on one more additional variable.

3.
We fix reference values for \(e_0\) and \(i_0\), and we introduce new variables \(\eta\) and \(\iota\) such that \(e=e_0+\eta\), \(i=i_0+\iota\).

4.
We expand the Hamiltonian in power series around \(\eta =0\), \(\iota =0\) up to order 3 in \(\eta\), \(\iota\).

5.
We split the resulting Hamiltonian into the linear part and a remainder. We compute the generating function and the canonical transformation of coordinates to remove the remainder to higher orders.

6.
Once obtained the new normal form, we disregard the remainder, so that the two actions corresponding to eccentricity and inclination become constants of motion.

7.
The initial values of the new constants of motion, which are the two additional proper elements, are obtained backtransforming the canonical transformations in terms of the original variables, namely in terms of the initial data.
For a specific case, we compute the osculating and proper elements by integrating the equations of motion and by computing the normal form using a Mathematica\(^{\copyright }\) program. We summarize below the steps of the procedure which will be implemented for each of the fragments of the case studies analyzed in the next sections.

Step 1. INPUT: set the normalization parameters: maxSteps=maximum normalization steps, maxR=number of terms kept in the remainder after each step, maxTaylor=maximum order of the Taylor expansion in the Lie Series, T=time span of propagation, step=integration step size.

Step 2. INPUT: Initialize the variables \((a, e, i, M, \omega , \Omega )\).

Step 3. Integrate Hamilton’s equations of the full Hamiltonian up to time T to get the osculating final elements.

Step 4. Compute the average of the Hamiltonian with respect to the mean anomalies of debris, Moon, Sun, and the sidereal time.

Step 5. Expand up to order 3 the averaged Hamiltonian (in the extended phase space) around the reference values \(e_0\), \(i_0\).

Step 6. Compute the generating function up to order maxSteps.

Step 7. Compute the new Hamiltonian using the generating function determined at Step 6.

Step 8. Compute the analytic solutions by determining the new coordinates as function of initial coordinates.

Step 9. Determine the two proper elements by integrating the analytic solutions over the given interval and dividing by the length of the interval.
Test cases: proper elements and data analysis
Let us consider a concrete case formed by, say, N fragments. In practical applications, our backtracing procedure is the following: (i) we take the (initial) orbital elements of all N fragments at time \(t=t_0\); (ii) we compute the initial proper elements from the initial orbital elements; (iii) we propagate all fragments up to a time \(t=T\) to compute the (final) osculating elements; (iv) through averaging and normal form, we compute the final proper elements from the final osculating data; (v) we compare the final osculating and final proper elements with the initial orbital and initial proper elements.
Since the proper elements are quasiintegrals of motion, we expect that they retain the main features both in the initial and the final phase, thus reconnecting much better to the original elements than the propagated osculating elements. Of course, the reconnection through the proper elements is more effective in those cases in which the final osculating elements get more dispersed over time, thus losing their link with the original data. Concerning step (v), beside making a visual inspection of the plots in the planes (a, i), (i, e) of (initial and final) osculating and proper elements, we apply data analysis techniques by using the KolmogorovSmirnov (KS) test and the Variance Equivalence (VE) test of the errors between the osculating and proper elements taken at the initial and final times. We also compute the Pearson correlation coefficients of initial vs. final osculating elements, and initial vs. final proper elements.
Such methods, borrowed from statistical data analysis, are briefly recalled as follows^{16}.

(S1) KolmogorovSmirnov test (KST) is a goodnessoffit test where the null hypothesis says that two datasets were taken from the same distribution, while the alternative hypothesis states that they are not taken from the same distribution. We used the predefined Mathematica\(^{\copyright }\) function KolmogorovSmirnovTest, which returns the pvalue of the statistical test. The pvalue has to be compared with a significance level \(\alpha\) (default is 0.05), null hypothesis being rejected for \(p < \alpha\).

(S2) Variance Equivalence test (VET) is a statistical tool that checks if the null hypothesis \(H_0\), that the variances of two data sets are equal, is accepted or not. Depending on the datasets’ distributions and the needed assumptions, one of the following tests is applied: “Brown Forsythe”, “Conover”, “Fisher Ratio”, “Levene”. We used a Mathematica\(^{\copyright }\) function called VarianceEquivalenceTest, that automatically chooses the most appropriate test and returns the pvalue and the conclusion of the test.

(S3) Pearson correlation coefficient, usually denoted by r, is used as a statistical measurement of the relationship between two onedimensional datasets. It is defined as \(r = \frac{Cov[X,Y]}{Var[X]Var[Y]}\) and gives a real number belonging to \([1,1]\), where 1 means a total positive linear relationship, 0 means no relationship, and \(1\) means a total negative linear relationship between the two datasets.

(S4) To visualize the data and to understand the main features of a distribution, one can plot the histogram of the dataset. This plot shows the frequency of each element from the set. This is a useful tool to compare the distributions of two or more data sets.
Among the cases available on the database “Stuff in Space” (http://stuffin.space/TLE.json), updating daily the orbit data from “SpaceTrack” (http://www.spacetrack.org/), we select the following test cases: Ariane 44lp, Atlas V Centaur, CZ3, Titan IIIc Transtage with a number of fragments equal, respectively, to 35, 164, 139, 41. We consider the following time intervals: 25, 50, 100, 150 years.
The outcome of the data analysis is summarized in Tables 2 and 3, where we provide the comparison between different elements. Table 2 gives the results, including the pvalues, about the KolmogorovSmirnov test and the Variance Equivalence test for the errors between osculating and proper elements at different final times. It is remarkable that both tests are always rejected, showing that the errors associated to the osculating and proper elements follow different distributions. Table 2 shows also the ratio of the root mean square errors of osculating versus proper elements, supporting that the errors associated to the osculating elements are larger than those associated to the proper elements. Table 3 gives the Pearson correlation coefficients of the initial and final, osculating and proper elements at different times.
In the supplementary material we detail the results for a fragment sample, that we take from Ariane 44lp; the supplementary material is aimed to help in reproducing the methods described in the present paper and, precisely, to compute the osculating elements at the initial and final times, to determine the normal form, to get the analytic solution and to compute the proper elements for the specific fragment. The same procedure can be implemented for the other fragments to get the results obtained in this work.
Using the data in Table 3, Fig. 2 summarizes the Pearson correlation coefficient between the initial data and the final osculating and proper inclination, as well as the initial and final proper inclination. In all sample cases, the correlation between the initial and final proper elements is always close to 1, while using the other sets we obtain discrepancies between the correlations of the initial and final states.
In the case of Titan IIIc Transtage there is a weak correlation between initial and final osculating elements, a better Pearson coefficient between initial osculating elements and final proper elements, and an almost perfect fit of initial and final proper elements. The other three sample cases have a similar behavior.
Atlas V Centaur
As an example, we analyze in detail the statistics of Atlas V Centaur and we describe the numbers shown in Tables 2 and 3. Table 2 shows that the KS test and the VE test are always rejected, both for eccentricity and inclination, at all times we investigated, namely 25, 50, 100, 150 years. Hence, the errors for osculating and proper elements follow different distributions with the errors associated to the osculating elements being larger than those of the proper elements.
The Pearson correlation coefficient in Table 3 tends to be constant when we compare the proper elements at different times. This result confirms the near constancy of the proper elements for a long period of time.
Figure 3 shows the evolution of the osculating elements in the plane ai compared with the evolution of the proper elements in the same plane (left); it also shows the distribution of the inclination (right) for the times 10, 25, 50, 100, 150 years in case of osculating (top) and proper (bottom) elements.
As it can be seen from the plots in Fig. 3, the osculating inclination starts to spread around 25 years; the spread increases with time. On the contrary, the proper inclination is kept almost constant, thus allowing to reconstruct the distribution at the initial time. This fact is also confirmed by the histograms and the associated Pearson correlation coefficients.
Titan IIIc Transtage
It is known that^{17} on February 21, 1992 an explosion of Titan IIIc Transtage produced several debris. All debris have been tracked and their coordinates at the present time can be found on “Space track”. We test our procedure, assuming to ignore the breakup time and propagating backward all fragments for a period of time equal to 29.5 years. The following results confirm the validity of the procedure based on the computation and comparison of the proper elements. In fact, like in the other cases, the KS and VE tests are rejected for all times with errors bigger for the osculating than for the proper elements. Beside, comparing the osculating elements at the present time and at the final time we obtain a Pearson correlation coefficient equal to 0.99886 for the eccentricity and 0.888539 for the inclination. On the other hand, comparing the proper elements at present and backward in time, we find a Pearson correlation coefficient equal to 0.999997 for the eccentricity and 0.999947 for the inclination.
Two mixed cases
We finally test our method by mixing the cases of CZ3 and Atlas V Centaur; the results are given in Fig. 4, which shows the evolution of the osculating and proper inclinations at different times (10, 25, 50, 100, 150 years). Through the proper elements, we succeed to distinguish two different clouds. In fact, while the osculating elements change so much that we cannot recognize the two groups after 50 years, the proper elements keep constant over all periods of time. This conclusion is supported by the comparison of the Pearson correlation coefficients and the histograms as in the right side of Fig. 4.
Conclusions and perspectives
The potentiality of the procedure of computing hierarchical proper elements for the space debris is witnessed by the four sample cases that we have analyzed in the previous sections and whose results are summarized in Tables 2, 3 and Fig. 2. Based on a solid mathematical method, our approach provides a reliable and effective way to connect the proper elements of a set of fragments to a specific breakup event. Even not knowing the breakup time, we can propagate backward the elements of the space debris to reconnect the debris to the parent body; we show the effectiveness of this procedure in the specific example of Titan IIIc Transtage. We have also shown a sample where the procedure allows one to distinguish between fragments generated by nearby breakup events. We stress that the computation of the proper elements can be extended to encompass limit cases of very small/large eccentricities and inclinations, as well as in the case of orbital elements close to resonances^{18,19}. Such cases can be dealt with a dedicated perturbation theory that generalizes our approach.
References
 1.
Klinkrad, H. Space Debris: Models and Risk Analysis (Springer, 2006).
 2.
Celletti, A., Pucacco, G. & Vartolomei, T. Proper elements for space debris. Submitted to Celestial Mechanics and Dynamical Astronomy (2021).
 3.
Hirayama, K. Groups of asteroids probably of common origin. Astron. J. 31, 185–188 (1918).
 4.
Knežević, Z. Asteroid Family Identification: History and State of the Art. In Chesley, S. R., Morbidelli, A., Jedicke, R. & Farnocchia, D. (eds.) Asteroids New Obs. New Models 318, 16–27 (2016).
 5.
Carruba, V., Aljbaae, S. & Domingos, R. Identification of asteroid groups in the z_{1} and z_{2} nonlinear secular resonances through genetic algorithms. Celestial Mech. Dyn. Astron. 133 (2021).
 6.
Knežević, Z. & Milani, A. Synthetic proper elements for outer main belt asteroids. Celestial Mech. Dyn. Astron. 78, 17–46 (2000).
 7.
Cook, G. E. LuniSolar perturbations of the orbit of an earth satellite. Geophys. J. 6, 271–291 (1962).
 8.
Gachet, F., Celletti, A., Pucacco, G. & Efthymiopoulos, C. Geostationary secular dynamics revisited: application to high areatomass ratio objects. Celestial Mech. Dyn. Astron. 128, 149–181 (2017) arXiv:1611.08916.
 9.
Rosengren, A., Amato, D., Bombardelli, C. & Jah, M. Resident space object proper orbital elements, in Conference paper, Maui, Hawaii, USA (2018).
 10.
Celletti, A., Galeş, C., Pucacco, G. & Rosengren, A. J. Analytical development of the lunisolar disturbing function and the critical inclination secular resonance. Celestial Mech. Dyn. Astron. 127, 259–283 (2017) arXiv:1511.03567.
 11.
Gkolias, I., Alessi, E. & Colombo, C. Dynamical taxonomy of the coupled solar radiation pressure and oblateness problem and analytical deorbiting configurations. Celestial Mech. Dyn. Astron. 132 (2020).
 12.
Kaula, W. M. Development of the lunar and solar disturbing functions for a close satellite. Astron. J. 67, 300 (1962).
 13.
Celletti, A. & Galeş, C. On the dynamics of space debris: 1:1 and 2:1 resonances. J. NonLinear Sci. 24, 1231–1262 (2014) arXiv:1408.1254.
 14.
Efthymiopoulos, C. Canonical perturbation theory; stability and diffusion in Hamiltonian systems: applications in dynamical astronomy. Workshop Ser. Asociacion Argentina de Astronomia 3, 3–146 (2011).
 15.
Poincaré, H. Les méthodes nouvelles de la mécanique céleste (GauthierVillars, 1892).
 16.
Cowan, G. Statistical Data Analysis (Clarendon Press, 1997).
 17.
Cowardin, H. et al. Observations of Titan IIIC Transtage Fragmentation Debris. In Ryan, S. (ed.) Advanced Maui Optical and Space Surveillance Technologies Conference, E5 (2013).
 18.
Lemaitre, A. & Morbidelli, A. Proper elements for highly inclined asteroidal orbits. Celestial Mech. Dyn. Astron. 60, 29–56 (1994).
 19.
Celletti, A., Gales, C. & Lhotka, C. Resonances in the Earth’s space environment. Commun. Nonlinear Sci. Numer. Simul. 84, 105185 (2020) arXiv:1912.04593.
Acknowledgements
The authors acknowledge EU H2020 MSCA ETN StardustReloaded Grant Agreement 813644. GP acknowledges the support of MIURPRIN 20178CJA2B “New Frontiers of Celestial Mechanics: theory and Applications” and GNFMINdAM. We thank the anonymous reviewers for valuable comments that helped to improve this work. A.C. and T.V. thank D. Marinucci for several interesting discussions on the statistical data analysis.
Author information
Affiliations
Contributions
A.C., G.P., T.V. contributed equally to this work.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare no competing interests.
Additional information
Publisher's note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Celletti, A., Pucacco, G. & Vartolomei, T. Reconnecting groups of space debris to their parent body through proper elements. Sci Rep 11, 22676 (2021). https://doi.org/10.1038/s4159802102010x
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s4159802102010x
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.