Re-entry predictions of space debris for collision avoidance with air traffic

The need for a safe and efficient integration of space vehicle operations into air traffic system to minimize the risk of impacts of spacecraft and aircraft and to sustain a steady air traffic is evident. This work provides a strategy toward more efficient management of uncontrolled re-entries by combining uncertainty propagation analysis with the FAA’s ConOps in the framework of the FAA’s Next Generation Air Transportation System (NextGen). The paper considers the scenario where a spacecraft re-entry is completely uncontrolled. During such re-entry, predictions are very difficult and are affected by various sources of uncertainty. Then the resulting position distribution throughout airspace boundaries is analysed and the impact on air traffic is estimated by defining protected airspace areas. The impact of the size of the protected area on the air traffic control is analysed and strategies for re-routing of air traffic are proposed. To verify the proposed results, the re-entry of the core stage of the CZ-5B R/B launcher is simulated.


Introduction
Since the beginning of the Space Age with the launch of Sputnik 1 in 1957, more than 24,400 catalogued orbiting objects have re-entered so far into Earth's atmosphere. Commonly, there are 200-400 objects crash into Earth annually, corresponding to about 50% of the total amount of returning mass. At altitudes between 75 and 85 km the re-entry vehicle (RV) may be broken into many pieces of debris because of aerodynamic and aerothermal loads. However this break-up depends on the materials that make up the object, so that re-entries of large intact objects can occur. There are two different types of re-entry : controlled and uncontrolled. During a controlled re-entry the target zone can be chosen to avoid any collision whereas during uncontrolled re-entry it is difficult to predict the re-entry timing and location [1]. Approximately, there are 70% of the re-entries of intact objects which are uncontrolled, corresponding to 50% of the returning mass. A data analysis of large intact objects, re-entered without control during 10 years from 2008 to 2017 shows that 4% of intact objects were upper stages with a mass exceeding 5 metric tons [2]. The most recent and significant case of uncontrolled re-entry of large intact object is the re-entry of the core stage of the rocket that launched Tianhe, CZ-5B R/B, in May 2021 [3]. This kind of atmospheric re-entry is more influenced by external disturbing forces, meaning that it is very difficult to predict their movements. Indeed various source of uncertainties, from phenomenons in the physical world or approximation models of the real world, affect the re-entry time and location predictions [4]. However, specific methods and procedures have been developed to provide a good assessment of re-entry trajectories. In fact, uncertainty propagation based on nonlinear Monte-Carlo (MC) simulations can be used in order to provide high-precision results [5,6]. On the other hand, the constant development of air traffic, the number of passengers carried and the lethality of each debris object for aircraft raise the need for a safe and efficient integration of RV 1 3 operations into the air traffic system [7,8]. Nevertheless, some methods used today based on restricted airspace contribute to inefficiencies for National Airspace (NAS) users, such as reroutes, delays, longer flight times leading to increase operating costs [9,10]. The main problem is that the size of the restricted area is often too large or the area is activated too early, resulting in false conflict detection. In order to avoid false conflict detection, the Federal Aviation Administration (FAA) published in May 2020 a Concept of Operations (ConOps) based on two different and new Airspace Management Methods to ensure Commercial Space Integration into the National Airspace (NAS) in the framework of the Next Generation Air Transportation System (NextGen) [11][12][13].
In this paper we propose a methodology for the integration of space traffic management (STM) and air traffic management (ATM) procedures. This methodology is based on the identification of the possible re-entry trajectories, affected by a certain uncertainty distribution, and their intersection with the air space occupied by air traffic. The size of the hazard area and the time at which the RV generates potential collisions are defined and, depending on the flight conditions of the airplanes within the NAS, appropriate diversions are proposed.
For the purpose of avoiding collisions between the RV and other NAS users, a hazard-area is defined at the entrance of the NAS, at flight level FL600. Then the time needed to clear this area is computed for each aircraft inside the area, and then it is compared to the time remaining before the RV entrance into the NAS. The purpose of this approach is to reduce the size of the restricted area, while minimizing the impact on air traffic. As long as the time to clear the hazard area is greater than the estimated time for the RV to reach the flight level under consideration, no diversion is needed for the air traffic. To achieve this, the dynamics of the atmospheric re-entry is computed and propagated from a given initial altitude up to the desired flight level [14,15]. The position dispersion obtained at FL600 is analysed in order to determine a polygon which contains all the possible impact points. The hazard-area to be evacuated is defined as the 95% error ellipse, enlarged with a buffer area in order to take into account aircraft separation standards [16][17][18][19].
The rest of the present paper is organized as follows: First, atmospheric model, Earth's shape and gravitational model are formulated briefly in Sect. 2. Next, the re-entry scenario, and uncertainty propagation are detailed in Sect. 3. These models are kept extremely simplified, to assess the overall procedure. In an operative scenario, more sophisticated models can be used to improve the accuracy of the re-entry prediction. Then, the airspace management methods are presented in Sect. 4. Finally, numerical results on one scenario obtained by the hybrid technique are given in Sect. 5 and the conclusion section ends this paper.

Problem statement
It is assumed that the re-entry vehicle (RV) during its trajectory is subject only to gravitational forces and aerodynamic drag. This section presents the atmospheric and gravitational models used in this work.

Atmospheric model
While gravitational forces are present throughout the whole trajectory, aerodynamic forces dominate at altitudes below 50 km. Furthermore, uncertainties in the aerodynamic loads can cause significant errors at impact. The atmospheric model must be able to represent the vertical distribution of pressure, density, and temperature. It is worth noticing that from pressure, density and temperature, there are many additional derived quantities such as dynamic viscosity, molecular mean free path length, and molecular collision frequency. In this paper, only the expression of density will be detailed. For more details, see [14]. As a starting point for developing a simple analytical model, let us assume that the altitude variable is the geopotential altitude h. By replacing the geometric altitude Z by the geopotential altitude h, we neglect gravity variations with altitude. Moreover, we assume that the atmosphere is entirely isotherm. Thus the density is expressed as: where 0 , T 0 , g 0 are the sea-level values of the air density, temperature and gravity, respectively. The gas constant, R, is set to 287 J/kg K. The lumped parameter RT 0 ∕g 0 also known as the atmosphere scale height H is 8.434 × 10 3 m. However, since the atmosphere is not isotherm through most of its altitude range the initial model has to be refined. As proposed by Regan and Anandakrishnan [14] we can adjust the two parameters 0 and H in order to obtain a good approximation of the U.S Standard Atmosphere 1976 for the altitude range of interest. For that purpose, an acceptable two-parameter atmosphere model can be written as follows: where and Equation (2) is very useful in obtaining closed-form solutions of re-entry vehicle trajectories. In Fig. 1 (4) H = 6.7 × 10 3 m between this model and the U.S. Standard Atmosphere 1976 is given. As expected, since a single reference height and reference density are used in this work, the model is reliable only until 120 km. Nevertheless the use of a more precise model for drag would improve the results still keeping the approach proposed in this article applicable.
According to the brief summary of some historical reentries given by [15], the re-entry point is chosen to be at an altitude of 80 km, so that this exponential atmospheric model fits well over the altitude range of interest, say from 5 to 80 km.

Earth's form and gravitational field
For our purposes, gravity can be defined as an interaction between a re-entry vehicle of mass m and Earth of mass M e , both treated as mass particles. Indeed, if we consider Earth as a sphere of uniform density, then the external gravitational field of Earth is the same as the acceleration obtained by concentrating all mass at the center of the sphere, i.e., a point mass. Thus Earth's gravitational acceleration can be expressed as: where G is the universal gravitational constant, R is the distance separating the Earth and the re-entry vehicle, R e and h are the radius of the Earth and the altitude of the RV above the surface of the Earth, respectively. For our altitudes of interest (from 80 km to the ground), the altitude h can easily be ignored in comparison to the Earth's radius R e [14]. As consequence, the Earth's gravitational acceleration will be constant, and set to 9.81 m/s 2 .

Scenario and initial conditions
The scenario considers an atmospheric re-entry of a space vehicle (SV) from near-orbital speed. Slightly before reentering the atmosphere, the spacecraft experiences dysfunctions of the control mechanisms so that the ability of controlling is lost. Therefore the re-entry becomes completely uncontrollable. The geometry of the RV is supposed to be cylindrical. The approximated dimensions and aerodynamics characteristics are based on the CZ-5B core stage, which re-entered Earth's atmosphere on May 8th 2021 after an uncontrolled re-entry [3]. These characteristics are summarized in Table 1.
Moreover, we assume that the RV is not capable of generating lift forces, so that it is only subjected to drag and gravity. As consequence, the restriction of planar motion results in no loss of generality, and can lead to closed-form expressions that can be very useful in assessing the reentry vehicle performances. The planar trajectory is shown in Fig. 2, where the following three sets of axis are used:   An inertial frame such that the X I Z I -plane contains the velocity vector V throughout the motion. 2. (X l , Y l , Z l ) : A local frame such that X l Y l -plane is that of the local horizontal and the Z l -axis is along the local vertical. 3. (X m , Y m , Z m ) : A moving frame attached to the RV such that the X m Z m -plane is coincident with the trajectory plane and the X m -axis is along the velocity vector at all times.

Re-entry vehicle particle mechanics
In this section, we present the equations of planar motion of the RV. We assume that the motion is determined by the initial conditions, the gravitational field and the atmospheric contact forces while luni-solar perturbations are neglected. Furthermore, we only consider a single stage, i.e, there is no break-up or explosion during the re-entry. The atmospheric flight phase of the re-entry trajectory can be described assuming a planar motion and no lift in a non-rotating environment according to the following equations: where is the ballistic coefficient, V the inertial velocity magnitude, the flight-path angle, h the altitude of the spacecraft and the downrange. The ballistic coefficient is defined as : where S is the reference area (usually the maximal cross-sectional area) and C D the drag coefficient. Figure 3 shows the typical altitude profile during the uncontrolled re-entry. This profile will be used in order to determine the time remaining before the RV reaches the National Airspace (NAS). In the following we consider that the atmospheric re-entry point is at an altitude of 80 km. Moreover, to simplify the procedure of describing the 2D position dispersion of the RV at a given airspace level, we define the X 1 -axis to be cross range (positive to the right from the direction of trajectory), the X 2 -axis down range, and the X 3 -axis up. Then we focus on the state vector X = (X 1 , X 2 , X 3 X 4 , X 5 , X 6 ) where X 1 , X 2 , are the positional coordinates and X 4 , X 5 , X 6 the velocity components. From Newton's Second Law, we can relate the three positional coordinates X 1 , X 2 , X 3 to the applied forces of aerodynamics and gravity. These equations will require six initial conditions : three positional and three velocity. The six equations of motion are given by: where is the density at the altitude X 3 and V the velocity magnitude which is given by:

Uncertainty propagation
Let us consider that the value of a typical trajectory is characterized by the proximity of its impact point to an

Fig. 3
Altitude profile during uncontrolled re-entry intended point (i.e, a target) such that, the closer the better. However, uncertainties in the initial conditions as well as in the physical properties of the atmosphere that Space Vehicle (SV) is reentering contribute to uncertainties in the impact point. A useful procedure is to define a nominal trajectory, which is assumed error free. Then, the uncertainty will be distributed around this nominal trajectory. In this paper we will only consider uncertainties in the initial conditions and in the drag coefficient. The variations of the initial state vector X 0 and the drag coefficient C D are performed by adding a randomized error to the respective nominal variables. A variety of physical sources can be responsible for these uncertainties. Even though it is not possible to assign a numerical values of arbitrary precision to these errors, it is possible to establish a mean value of a given error source and a distribution of likely values around this mean. In this paper we will consider that errors are given by a Gaussian distribution characterized by a standard deviation of likely values about the mean values (which are the nominal initial conditions and drag coefficient) [5]. Therefore each sample can be produced using a Gaussian pseudorandom number generator, such that errors on initial conditions and drag coefficient are considered as a zero-mean process about the nominal initial conditions and drag coefficient. The randomization of these variables is defined as follows: is the vector of the variables affected by uncertainties. Thus, each component is randomized with a standard normal distribution (mean equal to 0 and standard deviation equal to 1). These values are then multiplied with a standard deviation specific to each component. This uncertainty propagation is based on a standard deviation of 10 m for each position component, 10 m/s for each velocity component and 0.004 for drag coefficient, as suggested in the literature [14]. This Gaussian distribution is used to create a large number of sampled values [10]. Then the re-entry simulation is performed as a Monte-Carlo simulation in order to propagate these uncertainties, where the number of sampled values N samples is about 1000. The 1000 trajectories are determined by integrating the first ordinary differential equations Eqs. (11)(12)(13)(14)(15)(16). The 1000-trajectory Monte Carlo simulation are shown in Fig. 4. All trajectories are propagated until ground impact time.

95% ellipse error
As explained in Sect. 3.3, a variety of error sources contribute to discrepancy between the position of the intended target and the position of the RV at impact. Let the state vector associated with the nominal trajectory be designated by and the actual state vector by . Let us define the state error vector as: A precise statistical description of the error state is given by the state covariance matrix P which can be defined as the expectation of all possible pairs of error vector components: where P is a n × n matrix. For a particle representation, n is equal to 6. However for our interests we will restrict the error vector to 2D positional states only, in which case n is equal to 2. As consequence, P is a 2 × 2 matrix.
The error point position is the result of combination of coordinates uncertainties, thus the shape and size of the area that the point probably lies within its boundaries depends upon the nature of coordinates correlation and the value of uncertainties. Commonly, the position error takes an ellipsoidal shape with vary flattening according to the strength of the coordinates correlation [10,16]. Figure 5 shows the position dispersion for a 1000-trajectory Monte Carlo simulation. The simulation has been performed until the target altitude of 18km ( ≈ FL600).
It should be pointed out that the position error takes an ellipsoidal shape as expected. As consequence, it is Fig. 4 1000-trajectory Monte Carlo simulation relevant to define an error ellipse, which is an iso-contour, and a 2D confidence interval for a set of 2D normally distributed data samples (i.e initial conditions and drag coefficient distributions). For a given level of confidence, for example 95% , this ellipse defines the region that contains 95% of all impact points at FL600. Each error ellipse is defined by the following parameters: center position, heading and magnitude of semi-axes and level of confidence [17]. First, let us consider an axis-aligned ellipse, centered at the origin with a major axis of length 2a and a minor axis of length 2b. Its equation is defined as follows: where s defines the scale of the ellipse. This scale is associated with a given confidence level (e.g. a 95% confidence level corresponds to s = 5.991 ) [18]. Since our data are correlated (the covariance matrix P is not diagonal), the resulting error ellipse is not axis-aligned. This case corresponds to a rotation of the ellipse by an angle . As for the axisaligned case, the eigenvalues of the covariance matrix P represent the spread of the data whereas the eigenvectors represent the direction of the error ellipse axes. In other words, the eigenvalues represent the variance of the data in the direction of the eigenvectors [19]. The lengths of the semi-major axis and the semi-minor axis are given by: where min and max are the minimal and the maximal eigenvalues, respectively. Moreover, to compute the orientation of the ellipse, we calculate the angle between the largest eigenvector v max and the X-axis: where v x max and v y max are the components of v max along the X-axis and Y-axis, respectively. In addition, the ellipse is centered on the average X and Y positions of the data samples at FL600. Figure 6 shows the 95% error ellipse for the 1000-trajectory Monte Carlo simulation, at FL600(≃ 18 km) . The vectors shown by pink and green arrows in Fig. 6, are the eigenvectors of the covariance matrix of the data samples, whereas the length of the vectors corresponds to the eigenvalues. Now, we have to focus on the evolution of the 95% error ellipse parameters (i.e the lengths of the semi-major axis, a and the semi-minor axis, b) during the fall in order to characterize the area that has to be cleared to avoid collisions. For that purpose, the 95% confidence ellipse is computed for different target altitudes (i.e altitude at which the impact is considered). Fig. 7 shows the evolution of the semi-major axis and the semi-minor axis lengths, respectively, with respect to the altitude. The legends of both Fig. 7a and b give the time at which the ellipse has been determined, with origin of time at 80 km (i.e at t = 0 the RV is at an altitude of 80 km). It is worth noticing that graphics should be read in the direction of decreasing altitude.
The comparison of Fig. 7a and b shows that the semiminor axis is increased during the fall whereas the semimajor axis is decreased. Overall, this brings to an increase in the potential impact area (from 63 to 90 square kilometers). It is remarked that, since the trajectory becomes practically vertical below 40 km altitude (see figure 4), the initial error in the along-track component of the velocity

Airspace management methods
This section provides an overview of the Federal Aviation Administration's (FAA) methods for handling the growing needs of space operations in the framework of the FAA's NextGen program, which is the FAA-led modernization of America's air transportation system to make flying even safer, more efficient, and predictable [11].

Dynamic airspace separation
As Launch Vehicle/Re-entry Vehicle (LV/RV) operations are expected to increase across the National Airspace System (NAS) the FAA developed methods for ensuring a cost-efficient integration of SV into the NAS. Currently, large amounts of airspace are restricted for every LV/RV operation and for a long time, which increases costs for NAS users. To reduce the impact of LV/RV operations, the FAA's NextGen office recently proposed two more efficient separation concepts for LV/RV operations called Space Transition Corridors (STC, a method of airspace segregation) and Four-Dimensional Trajectory Deconfliction (a method of separation) [12]. These new concepts are based on the definition of protected volume which encapsulates potential hazards associated with the RV location and moves. This aggregate volume is derived from two other volumes: • The first volume protects against collisions and encounters between the RV and other NAS users, using the vehicle position uncertainty of the nominal operation. • The second volume is a separate volume used to protect aircraft from a hazard in the immediate aftermath of an off-nominal event (e.g., breakup, explosion, or thrusttermination).
As we focus on a single stage, with neither break-up nor explosion, we only consider the first volume, which mitigates the vehicle position uncertainty. This restricted airspace volume which covers a vertical area from the surface to an unlimited altitude (e.g FL600), must remain clear over the duration of the reentry operational window. On the other hand, the current approach for determining the volume of airspace that has to be segregated focuses on mitigating the danger of falling debris to aircraft. This concept is based on calculating the probability of falling debris causing a casualty on an aircraft [8]. This approach is sufficient for the definition of Temporary Flight Restrictions (TFR). However, as Space Transition Corridors and Four-Dimensional Trajectory Deconfliction call for more dynamic airspace separation, other metrics have to be taken into account such as the time to clear the hazard volume, controller workload, possibility of loss of radar reparation between aircraft during stressful situations, and probability of Near Mid-Air Collision (NMAC). In the following, our scenarios are focused on metrics related to the concept of Space Transition Corridors. Furthermore, these STC can be dynamically activated and deactivated on a just-in-time basis to enable the RV to transition safely through NAS while minimizing the effect on other NAS operations [12]. For the purpose to implement this dynamic activation, we need to determine the time remaining before the RV reaches the NAS as well as the time needed to clear the restricted area at a given airspace level (e.g FL600), so called the look-ahead time. Therefore, based on the air traffic at FL600, a strategy to clear the hazard-area for this specific altitude will be defined.

Look-ahead time to clear the hazard-area
As consequence of the huge number of dynamically changing variables, the need for an automated decision-support system to assist Air Traffic Controllers (ATC) in determining the protected area with enough lead time to clear atrisk aircraft is evident. This lead time is highly dependent on the scenario, the air traffic density, the type of aircraft involved and the size of the protected area. Thus, to implement these new separation concepts, based on dynamic protected volumes, accurately predicting a parameter called the NAS Automation Boundary Entry Time (NABET) is required in order to avoid false conflict detection alerts. The NABET is the time at which the vehicle will enter the NAS automation boundary (i.e FL600 in our case) with a high standard of accuracy [9,12]. However, due to uncertainty propagation the time at which the RV will re-enter the atmosphere is not exactly known. Therefore, an uncertainty window associated with a confidence level should be defined. In view of the statistical distribution of the predictions, an uncertainty window amplitude of ±30% around the estimated nominal re-entry time guarantee a confidence level ≥ 95% [4]. Thus, before the NABET is received, re-entry time uncertainty may result in false conflict alerts, causing maneuvers based on false alert and then additional costs for NAS users. Indeed, during the descent ATC must maintain minimum separation between aircraft and the RV as well as aircraft and aircraft. In our scenarios, ATC is assumed to know the aircraft flight profile beforehand. Moreover, there are automation tools that can provide ATC with knowledge of when an expected conflict would occur. Then, once the ATC is aware of a conflict, other NAS users are maneuvered as necessary to resolve conflicts involving the RV, as the RV is not capable of complying with ATC tactical control instructions (uncontrolled re-entry). Once the NABET is computed, the associated envelope of protected airspace can be accurately determined in order to avoid collisions. Then the protected airspace is activated, so that ATC have to direct aircraft that are in the protected area to exit as quickly as possible and aircraft outside the protected volume to avoid entering. Thus, the NABET lead time (i.e the time at which the NABET is received) needs to provide enough time for conflicts in the deconfliction horizon to be safely resolved before the RV enters the airspace (i.e FL600). It should be emphasized that both the NABET and the NABET lead time are unique to a given operation and influenced by multiple factors (e.g., location of operation, vehicle speed, etc). ATC and pilot response times have also to be taken into account for computing the time needed to clear the restricted area. This response time depends on the ATC workload, so on the number of aircraft that are within the ATC sector.

Hazard-area evacuation
This section provides the different concepts used by ATC in our scenarios to clear the hazard-area : the geometric form of the restricted area as well as the proposed evacuation process.
In general, hazard-areas involving LV or RV are modeled as convex polygons. In our study, the 95% error ellipse is taken as the original hazard-area. Then the algorithm calculates the axis lengths, of an enlarged ellipse encompassing the original ellipse with buffer areas added to take into account separation standards. However, today because of the lack of separation standards, procedures and/or techniques available to controllers to separate aircraft from LV and RV, there is no true separation standards similar to today's standards for civil aircraft. Further studies are needed to develop generalized spacecraft categories regarding design layouts, flight performance and operational profiles [7]. Thus, in this work the horizontal separation distance used for the RV is 5 nautical miles (NM) (IFR separation with radar). In the re-entry scenario that will be discussed the vehicle experiences dysfunctions of the control mechanisms. Therefore, the re-entry becomes completely uncontrolled which requires a restricted area to be activated. At this time, the controllers' radar displays show how much time the controllers have left to get aircraft out of the newlycalculated area to avoid potential collision with the RV. For that purpose, hazard-area will be studied at a specific altitude to compare the time needed to clear this area and the time remaining before aircraft impact. This specific altitude is chosen to be the top ceiling of the NAS (i.e FL600). The aircraft traffic used in these scenarios are based on a uniform random distribution of aircraft inside the hazard-area at this specific altitude. Given the size of the 95% ellipse and the 5 NM (i.e.: 9.26 km) separation standards, there will be only three aircraft inside the hazard-area at the beginning of the simulation. Then ATC instructions are given to aircraft to remove flights from the hazard-area. Probabilistic times for ATC and pilots to issue, respond, and execute commands (e.g., due to pilot read-backs and delay until maneuver initiation) are used, and are based on prior MITRE Corporation research on controller workload [6]. The challenges of these scenarios are to compute both times to clear the enlarged hazard-area and time before impact, in order to provide separation between aircraft and the STC at the entrance of the NAS. The hazard-area evacuation algorithm calculates the required change in the aircraft heading to exit from the ellipse as quickly as possible. Under the assumption that aircraft are moving inside the volume at a constant speed, in a straight line after finishing the turn, the time needed by the aircraft to exit the ellipse can be computed. Indeed, for aircraft with a true airspeed greater than 170 kts, the rate of turn is given by: where g is the gravitational constant, v the aircraft true airspeed and the bank angle [6]. Furthermore, the time duration T of the heading change is: Thus, the total amount of time it takes for each aircraft to exit the ellipse can be calculated as the sum of the time duration of the heading change and the time in straight line after finishing the turn, so that the turn angle corresponding to the shortest total amount of time will be used as the command angle of turn to the aircraft to exit the hazard-area as quickly as possible [6]. However, in our scenarios, there are multiple aircraft in the polygon, so ATC needs to command aircraft to exit the polygon as quickly as possible while avoiding conflicts between these aircraft. For that purpose, the aircraft that takes the least amount of time to exit the polygon is commanded to execute the heading change calculated by the algorithm. Then if a conflict occurs with the two remaining aircraft, these aircraft are maneuvered to avoid the conflict with the first aircraft. Then, if a conflict between these two aircraft appears, the algorithm will solve the conflict by maneuvering only one of these two aircraft. To achieve this goal, for each aircraft the algorithm computes from the list of possible heading changes (i.e heading changes that don't generate a new conflict with the first aircraft) the heading change that allow to avoid the conflict while minimizing the delay caused by the maneuver. Then the shortest time it takes to exit the hazard-area is selected, the corresponding aircraft is maneuvered with the associated heading change. If no solution is found, each aircraft continues on its nominal trajectory.

Time to clear the hazard-area
As explained in Sect. 4.2 current methods of calculating the restricted area may result in false conflict detections due to the large amounts of airspace that are restricted and for a long duration. In this paper, we propose a technique to find the altitude from which the time required to clear the airspace will be greater than the time remaining before impact. The proposed approach can be summarized as follows: 1. Initialisation: Choose the altitude decreasing step dh, the uncertainty distribution, the initial velocity. The initial altitude is set to 80 km. 2. Uncertainty propagation: Propagate the trajectory until FL600 using Monte-Carlo simulations. Compute the time remaining before impact using the 1000-trajectory Monte Carlo simulation, while computing the time to clear the hazard-area using the air traffic configuration and the hazard-area evacuation algorithm presented in Sect. 4.3. 3. Comparison: Compare both times. If the time remaining before impact is greater than the time to clear the hazard-area, then the initial altitude is decreased by dh, the velocity reached by the RV at this altitude during the nominal trajectory is used as initial velocity for the next iteration. The main purpose of this approach is to reduce the impact on air traffic during the re-entry. Indeed, as noted by Regan and Anandakrishnan [14], the error ellipse grows in size throughout the re-entry segment, indicating that uncertainties attributable to the environment will make further contributions to positional errors, and that velocity errors translate into positional errors over time. Consequently, the process described above aims at reducing the size of the hazard-area at the entrance of the NAS. A flow diagram, illustrating the proposed algorithm is shown in Fig. 8.

Results and discussion
To illustrate the ideas and techniques proposed above, numerical results on one possible scenario are presented in this section. In this example, the air traffic is represented by random positions inside the 95% error ellipse as well as by random headings for each aircraft. First the hazard-area evacuation algorithm computes for each aircraft the heading corresponding to the shortest time to exit the area. Then, based on the time needed to clear the hazard-area, an iterative algorithm is applied to find the time when the hazardarea has to be activated to avoid collision between the RV and other NAS users.

Hazard-area evacuation scenarios
First, it should be emphasized that the air traffic generated by the algorithm is deconflicted (i.e 5 NM of horizontal separation) at the beginning of the simulation. Furthermore, given the lengths of the error ellipse's axis, there are only three aircraft inside the 95% error ellipse at the beginning of the simulation. In this example, the traffic will be studied at the specific flight level FL600. The aircraft true airspeed is set to 448 kts which corresponds to the average cruise speed of an A320. Furthermore the bank angle is set to 67 • (i.e the maximal bank angle for an A320 under the Normal Law). For each aircraft, the turn angle is considered to be in −60 • , 60 • . Both semi-major and semi-minor axes of the 95% error ellipse have been enlarged by 5 NM to take into account separation standards between aircraft and the RV. One may note that it is this enlarged ellipse that has to be cleared.
In this scenario, ATC's instructions to aircraft don't lead to any conflict so that each aircraft can follow ATC's instructions in order to exit the volume as quickly as possible. Figure 9 shows aircraft operating in the area at the beginning of the simulation for this scenario. Figure 9a only shows nominal trajectories and exit points before ATC's instructions, whereas Fig. 9b shows the new headings and exit points after ATC's instructions. The red vectors and points represent the velocity vectors and exit points before ATC's instructions, respectively, whereas the cyan vectors and points represent the velocity vectors and exit points after ATC's instructions, respectively.  The angles and rotations shown in Fig. 9 are measured in the horizontal plane and counter-clockwise. Table 2 summarizes the relevant metrics for this scenario.
Since these maneuvers don't lead to any conflict, each aircraft can follow ATC's instructions to clear the volume as quickly as possible. Consequently, in this scenario, the hazard-area is cleared 46.5 s after ATC's instructions whereas without ATC's instructions the total amount of time to clear the hazard-area is 51.7 s. Note that times computed here don't take into account ATC and pilot response times. Nevertheless, it is important to point out that the RV is considerably higher (in altitude) than all aircraft at the beginning of the simulation, and would not reach aircraft altitude for a few minutes (time to descent from 80 to 18 km ≈ FL600 ). This provides time for ATC to issue instructions to flights to clear the area before the RV reaches the NAS.

Earth's atmospheric re-entry
An Earth atmospheric re-entry similar to the CZ-5B reentry [3] is analysed, in particular the time remaining before impact, to find the RV's altitude at which the hazard-area has to be cleared to avoid collision with air traffic. The approach used is given by Fig. 8. For each iteration the time needed to clear the hazard-area has to be determined. As mentioned in Sect. 4.3 times for ATC and pilots to issue, respond, and execute commands are needed. These times are taken from prior MITRE Corporation research on controller workload [6]. Based on a probabilistic approach, they found that the time until the first reroute was given to the pilot was about 25.09 s for an average number of four aircraft inside the hazard-volume. Given our similar traffic density, results from [6], a fixed delay of 30 s is used to represent the time needed to the pilot to respond to ATC's instructions. This delay is added to the time needed to execute the maneuvers presented in Sect. 4.3. Therefore, for each aircraft the total amount of time (i.e pilots and ATC response time added to the time to execute maneuvers) is the time needed to clear the hazard-area.
The first interest of this approach is to focus on the evolution of the size of the restricted area with respect to the initial altitude. Figure 10 describes the evolution of the surface of both 95% error ellipse and the enlarged ellipse with respect to the initial altitude.
As shown by Fig. 10 when the initial altitude z 0 decreases, the surface of the hazard-area decreases as well. As consequence, at each iteration of the algorithm presented in Fig. 8, the initial altitude of the Monte Carlo simulation is decreased and consequently the surface of the error ellipse.  This procedure allows to reduce the number of aircraft involved in a potential conflict with the RV.
To achieve this procedure, we need to know at which altitude the time to clear the hazard-area is greater than the time before impact. When the RV reaches this altitude (that we call decision altitude), the corresponding hazard-area has to be activated and aircraft maneuvers have to begin to avoid collisions. The evolutions of both time to clear the hazardarea t c and time before impact t i with respect to the initial altitude z 0 are depicted in Fig. 11. It is worth noticing that the time needed to clear the hazard-area (blue line) seems to be constant from a given altitude. This is because from this altitude the size of the 95% error ellipse is negligible compared to the 5 Nm buffer so that the dimensions of the enlarged ellipse are mainly dominated by this buffer. The intersection between the red and blue lines gives the decision altitude, which is at 22 km in our simulation.

Conclusion
A new iterative procedure is being proposed to minimize the impact of uncontrolled re-entries on the NAS. Thereby, scenarios of Earth atmospheric re-entries have been analyzed to present separations standard used in this work. The 95% error ellipse, defined through the dispersion of the 1000-trajectory Monte Carlo simulation at the entrance of the NAS and the covariance matrix, is enlarged using these separation standards to compute the hazard-area that has to be cleared. Then the proposed algorithm computes and compares both time to clear the hazard-area and the time before impact. The procedure is performed until the RV reaches entrance of the national airspace, at FL600 ( ≈ 18 km). Given our simulations, the hazard-area has to be activated when the RV reaches an altitude of 22 km. This work is a preliminary step towards further interactions between uncontrolled re-entry and air traffic management optimization.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. No funds, grants, or other support was received.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.