Reconnecting groups of space debris to their parent body through proper elements

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 break-up. In this work, we present a procedure to back-trace 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 break-up, 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, CZ-3, 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.

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 break-up 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 break-up 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 ω and longitude of the ascending node . 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 www.nature.com/scientificreports/ body. Knežević and Milani 6 introduced also the synthetic proper elements based upon a numerical integration, a digital filtering of the short-period 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 break-up 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 break-up events of the satellites Ariane 44lp, Atlas V Centaur, CZ-3, 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 back-propagate the debris for a given time and compare the osculating or proper elements at the initial time and at the back-propagated 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, Kolmogorov-Smirnov 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 non-spherical 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 break-up 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 area-to-mass ratio lower than 0.74; (ii) there does not exist a public catalogue that provides information about the area-to-mass 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 non-spherical 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). 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: 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 time-scales, 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, semi-fast 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 (I, ϕ) are action-angle variables with (I, ϕ) ∈ B × T n , where B ⊂ R n is an open set and n denotes the number of degrees of freedom. In (2) the function H 0 (I) is the integrable part, ε ∈ R is a small parameter, H 1 (I, ϕ) 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 ε 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 H 1 can be expanded in Fourier series as where K ⊆ Z n and b k are functions with real coefficients. Let χ be the generating function of the canonical transformation from the variables (I, ϕ) to the new variables (I ′ , ϕ ′ ) given by where the action of S ε χ is defined by with {·, ·} the Poisson bracket operator. We determine S ε χ by requiring that the new Hamiltonian and H 2 is the remainder term of order ε 2 . Inserting the change of coordinates in (2), one obtains the transformed Hamiltonian which takes the desired form (3) provided χ satisfies the following normal form equation: Expanding χ in Fourier series, denoting the frequency by ω 0 = ∂H 0 ∂I ′ , one obtains that the generating function is given by the following formula, which is valid under the non-resonance assumption k · ω 0 � = 0: A higher order normal form is obtained by iterating the above procedure.
Recalling that the space debris model described above depends on fast, semi-fast 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 semi-fast (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 quasi-integral of motion for the averaged approximated model. After averaging over the mean anomalies and the sidereal time, we end-up with a Hamiltonian function with three degrees-of-freedom 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 . www.nature.com/scientificreports/ a canonical change of variables through a Lie series normalization, implemented through a Mathematica © 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 back-transforming 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 M depends on time, the Hamiltonian resulting from step 1 depends on (e, i, ω, �, 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 η and ι such that e = e 0 + η , i = i 0 + ι. 4. We expand the Hamiltonian in power series around η = 0 , ι = 0 up to order 3 in η , ι. 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 back-transforming 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 © 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 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 back-tracing 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 quasi-integrals 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 Kolmogorov-Smirnov (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) Kolmogorov-Smirnov test (KST) is a goodness-of-fit 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 © function KolmogorovSmirnovTest, which returns the p-value of the statistical test. The p-value has to be compared with a significance level α (default is 0.05), null hypothesis being rejected for p < α. www.nature.com/scientificreports/ assumptions, one of the following tests is applied: "Brown Forsythe", "Conover", "Fisher Ratio", "Levene". We used a Mathematica © function called VarianceEquivalenceTest, that automatically chooses the most appropriate test and returns the p-value 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 one-dimensional datasets. It is defined as 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.
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 p-values, about the Kolmogorov-Smirnov 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, Table 2. The p-value of the Kolmogorov-Smirnov test and Variance Equivalence test for the errors between osculating (eoe) and proper elements (epe) at different final times (25, 50, 100, and 150 years). The last two columns contain the ratio of the root mean square (RMS) between the osculating errors and the proper elements errors. The last row refers to an example where we back-propagate in time. A 0 p-value means a number lower than 10 −40 . www.nature.com/scientificreports/ 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. Table 3. The Pearson correlation coefficient obtained propagating forward in time between the initial osculating elements (ioe) and the final osculating elements (foe), between the initial osculating elements (ioe) and the final proper elements (fpe), between the initial proper elements (ipe) and the final osculating elements (foe), and between the initial proper elements (ipe) and the final proper elements (fpe) for eccentricity (ecc.) and inclination (incl.) at different final times (25, 50, 100, and 150 years). The last row refers to an example where we back-propagate in time.  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 a-i 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 break-up 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 CZ-3 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