Computation of Artificial Meteors Trajectory and Ablation

This work has been motivated by the lack of meteor data, which hinders the determination of the mass and composition of natural meteors. To advance the knowledge of meteors science, the Japanese start-up ALE Co. Ltd. (standing for Astro Live Experiences) designed hundreds of 1-cm diameter particles composed of materials representative of natural meteors, which will be release at an altitude of 300 km by an original payload in 2023. During entry, the light emitted by the particles will be analyzed by spectrometers and cameras on the ground to characterize the thermal and optical properties of the materials. This paper carries out a sensitivity analysis and aims to discuss the modeling parameters predicting full disintegration of the particles before touching ground. This paper presents the efforts undertaken to compute the trajectory of artificial meteors during their entry into Earth’s atmosphere. Various materials representative of natural meteors and their thermal response were modeled. The trajectory of artificial meteors was computed by solving the equation of motion including mass conservation. The influence of the drag coefficient, heat transfer coefficient, and geometric parameters was evaluated through the statistical analysis of Grid Sampling and Monte Carlo simulations. They were shown to have a sizeable effect on the trajectory, mass loss, latitude and altitude of demise. The computations demonstrate full demisability of the artificial meteor above 60 km, which is compliant with international safety regulations.


Introduction
Natural meteors science remains an active field of research as it enables the understanding of the origin of life. Of paramount importance is the meteor brightness efficiency, which relates the observed luminosity to the rate of variation of its kinetic energy (including velocity decrease and mass loss) and has been the subject of significant research over the last decades [1][2][3]. Meteors bright tails are also considered by the community as natural indicators of the phenomena occurring in the mesosphere [4], which is a layer whose properties remain unknown while affecting the upper and lower layers. The analysis of the tail signature provides valuable information such as the surrounding temperature, the electron concentration, and the composition of the materials. The information retrieval from these analyses is not straightforward and requires iterative algorithms coupling complicated physico-chemical phenomena such as emission, ionization, ablation, etc. A result of the complexity of the analysis is the significant uncertainty lying within the brightness efficiency [5,6]. To tackle this issue and advance meteor science, ALE Co. Ltd. (which stands for Astro Live Experiences) develops with its scientific partners a technology aiming at accurately releasing hundreds of meteors particles of known composition, mass and velocity. The analysis of the particle signature, complemented with high-fidelity simulations, representative ground experiments and comprehensive observation campaigns will provide calibrated, accurate and on-demand data to the community. The paper objectives are fivefold: (i) model the reentry mechanisms, (ii) implement the governing equations into a reliable simulator, (iii) compute the trajectory and mass loss, (iv) carry out a sensitivity analysis on the aerodynamic and mechanical parameters, (v) conclude on the safety of the mission, i.e. ensure the meteor completely vanishes after emitting light. The chapter is organized into four sections. Section 2 briefly introduces the proprietary technology developed by ALE Co. Ltd. to generate artificial meteors. Section 3 describes the governing trajectory and conservation equations as well as the physico-chemical and numerical modeling relevant for Earth's atmosphere entry. Section 4 describes the orbit propagation algorithms developed throughout this research. Section 5 carries out a sensitivity analysis of the parameters involved through a Monte Carlo analysis for a case scenario formulated by ALE Co. Ltd. and Sect. 6 discusses the numerical results. Finally, Sect. 7 presents the conclusions obtained.

Artificial Meteors
Meteor data are scarce, which hinders the accurate determination of space debris demise altitude and the advancement of meteor sciences, specifically in the field of material ablation and emission. For instance, the heat shield systems of (re)entry spacecraft and rockets are equipped with materials which encounter ablation when subject to high heat fluxes. Understanding ablation in meteor sciences will thus 1 3 The Journal of the Astronautical Sciences (2022) 69:  advance the knowledge of heat shield ablation, determine more accurately their size and purchase cost.
Determining the mass, composition and entry velocity of natural meteors from observation is complicated and warrants iterative procedures. To tackle this issue, ALE developed a unique shooting star technology based on the release of hundreds of particles of known mass, composition, trajectory from a constellation of satellites. Figure 1 illustrates the meteor-generation principle pioneered by ALE Co. Ltd. A satellite, depicted in Fig. 2, carrying hundreds of 1 cm diameter spherical meteor sources is launched. Upon altitude reached, the meteor sources are emitted at a specified position, direction, and velocity and become artificial meteors because of the friction of the atmosphere layers. The initial launch orbit of the ALE-1 satellite is an elliptical orbit with a perigee altitude of 485 km and an apogee altitude of 492 km. Then, the satellite will decrease their orbit thanks to its propulsion systems  To enable the reentry of the meteor source into the Earth's atmosphere, it is necessary to release the meteor source in the direction opposite to the satellite's traveling direction. The release system performance was assessed in vacuum chamber. The brightness of the meteor was also characterizing by testing the material in plasma wind tunnels operated at representative heating conditions [7]. ALE ondemand meteor data serve as a benchmark to assess the performances of the simulation tools, calibrate the cameras on ground and further characterize the complicated phenomena governing meteor brightness (emission, non-equilibrium gas chemistry, ablation, etc.).

Governing Equations
A common approach used in the community is to solve the equations of motion of the sphere from release from the satellite to its full ablation during Earth's entry. The model developed by Tokyo Metropolitan University [8] used a fourth order Runge-Kutta algorithm to solve the equations of motion as well as the mass loss equation as derived in Eqs. (1.a) and (1.b).
Equation (1.a) describes the linear momentum equation, where is the position vector, atm = − ( × ) is the velocity relative to a co-rotating atmosphere, with the inertial velocity vector and the angular velocity vector of the Earth, m is the mass of the meteor source, C d the drag coefficient, S the cross-section area of the meteor, ρ the atmospheric density and μ the Earth's gravitational parameter. Equation (1.b) corresponds to the mass loss equation due to ablation such that L * = 10 6 J ⋅ Kg −1 is the ablation specific heat and C h the heat transfer coefficient. Note that the linear momentum equation has been simplified by assuming a uniform mass transfer around the surface of the meteor [9]. Indeed, the rate of change of momentum turns out to be where l is the velocity of the material leaving the meteor with respect to the meteor. Considering a uniform mass transfer, this velocity becomes equal to zero. From arc-jets experiments [7] and modeling from the literature, ablation is not spherical. While a material response modeling of the particle ablation is beyond the scope of this work, the shape change was modeled with a cross-section area, as derived in The Journal of the Astronautical Sciences (2022) 69:1347-1374 the following section. This code yields the trajectory, the heating rate as well as the meteor source speed decrease due to friction (modeled by the coefficient C d ) as well as mass loss and radius decrease due to ablation (modeled by the coefficient C h ), as displayed in Fig. 3.

Cross-Section Area
The aerodynamic area changes as the mass decreases. Their relation is expressed by the shape change coefficient ν [10] as where S 0 and m 0 are the initial values before ablation of the cross-section area and mass respectively. The shape change coefficient ν is a parameter that characterizes the rotation of the meteor along the flight. In this study we consider ν = 0.66, that is, a meteor remaining spherical.

Atmospheric Parameters
A body moving close to the Earth is affected by the atmospheric drag. Satellites in Low Earth Orbit (LEO) encounter atmospheric drag from gases in the thermosphere (~ 80 km to 500 km), while objects like meteors are specially affected during reentry (~ 120 km). This perturbing acceleration depends on the atmospheric density ρ air and atmospheric temperature T air which need to be modeled with high accuracy in order to properly describe the constant change of the environment. For this study, the NRLMSIS-00 empirical model of the atmosphere [11] has been selected. Note that an updated model NRLMSIS 2.0 was recently developed [12], which concludes that the assessment of the new NRLMSIS 2.0 model does not seem to provide any significant improvement in comparison to its predecessor NRLMSIS-00. The NRLMSIS-00 enables us to use a single model by covering the atmosphere from the surface to lower exosphere (~ 0 km to 1000 km). The model ingests the F 10.7 and A p indices and outputs the total mass density, the temperature and the oxygen number density. The greater the Meteor source a speed, b heating and c radius decrease during entry 1 3 solar activity and geomagnetic activity are, the larger the atmospheric density and temperature of the air are. In this case, we select F 10.7 = 150 and Ap = 4, which corresponds to a medium activity.

Surface Temperature
The temperature of the surface of the meteor is a parameter that is taken into account in the atmospheric resistance, in particular in the drag coefficient. The meteors are acted upon not only by the aerodynamic heating, but by complex processes as melting and spalling. Consequently, the surface temperature becomes difficult to estimate. One option, which is selected in this article, is to consider it equal to the atmospheric air at the corresponding altitude:

Drag Coefficient
It is the main parameter in the atmospheric drag, which in the case of reentry objects is not constant, but significantly fluctuates along its trajectory. The expression considered for this parameter is given by [13]: such that: is the Mach number with R gas = 8.3144598 J ⋅ mol −1 ⋅ K −1 the molar gas constant, W air = 0.0289644 kg ⋅ mol −1 the molar mass of dry air, and = 1.4 the heat capacity ratio of air.
is the Reynolds number with R m the radius of the meteor and μ air the atmospheric viscosity given by [14]: The Journal of the Astronautical Sciences (2022) 69:1347-1374 with T 0 = 273.11 K , ref = 6.7894 ⋅ 10 −5 kg ⋅ m −1 ⋅ s −1 and T ref = 110.56 K reference values for viscosity and temperature respectively.

Heat Transfer Coefficient
It is the main parameter in the mass loss equation and it is related with the mass reduction of the meteor. The formula used in this study is taken from [15], [16,Chapter 5]: Equation (9) provides a formal definition of C h , although the integrands can be solved analytically, thus providing a more convenient expression.

Visual Magnitude
This is a proxy for the luminosity of the meteor during its atmospheric reentry. Due to the heating produced by the atmospheric drag during the reentry, the meteor emits energy in different regimes of the electromagnetic spectrum, including the visual range. Therefore, for the proposed experiment it is interesting to be able to quantify the expected luminosity of the meteor during the atmospheric flight, which is done adapting the approach presented by Dias et al. [17]. The traditional luminosity equation for a non-decelerating body (a reasonable hypothesis, since the negligeable residence time of the fluid around the asteroid, as compared with its orbital motion, allows for a quasi-stationary analysis of the flow) relates the meteor luminosity I to the dissipated amount of kinetic energy, assuming the latter is translated to the meteor surface as heat: where v ∞ represents the velocity relative to the atmosphere and is the unitless luminous efficiency for the considered bandwidth of the spectrum, which represents the amount of kinetic energy transformed into radiation in a specific bandwidth. For the visual bandwidth, the value = 10 −3 can be assumed. Once the meteor luminosity is known, the actual radiative spectral flux observed from the meteor by an observer located on ground can be estimated as: where d is the distance from the observer to meteor (we assume it coincides with the altitude over the Earth surface). Note that this definition of distance is different from astronomical observations where distances are measured relative to a fixed ground observer, where the station is placed. Right now the mission does not have any prefixed ground station and consequently we consider the distance to the Earth's surface as the natural distance available to apply the equation. The reader must take into account that the visual magnitude (M v ) defined here is related to an observer located at the subsatellite point. Finally, the luminous magnitude observed in the visual bandwidth and measured in mag0 units can be obtained from the following relation:

Initial Conditions
A priori initial conditions for the satellite and the meteor need to be introduced to properly propagate the orbit of the meteor. The present study is an illustrative example of the meteor's trajectory. Consequently, representative values for the problem are chosen for the sake of simplicity: instead of considering a sun-synchronous orbit, a satellite with a polar circular orbit and an altitude with respect to the equator of 375 km is considered. The meteor is shot just behind the satellite with relative velocity of 350 ms −1 . The release point of the meteor is defined by the orbital position of the mothership satellite, described through its orbital inclination (i), semimajor axis (a), eccentricity (e), longitude of the ascending node (Ω), and argument of latitude (u), which is the sum of the argument of periapsis (ω) and the true longitude (θ). These orbital elements for the satellite's orbit are: The start date is 2020/01/01, 00:00:00 UTC such that the orbit is propagated until one of the following conditions is satisfied: • The mass reaches a value equal or smaller than 10 -7 kg. • The meteor reaches the surface of the Earth.
The initial position of the meteor coincides with that of the mothership satellite from where it is to be released, and only the orbital velocity varies. The meteor is shot from the mothership with a relative velocity vector rel such that the meteor is released in orbit with a lower orbital velocity in order to ensure a reentry trajectory. This relative velocity vector is defined by its magnitude v rel and the yaw ( ) and pitch ( ) angles as defined by the Tait-Bryan rotation sequence (Z-Y-X) with respect to an orbital frame defined with the X-axis along the orbital velocity vector of the mothership satellite, the Z-axis along the nadir direction, and the Y-axis completing a right-handed frame. With respect to this orbital frame, the release conditions of the meteor relative to the mothership satellite are such that the meteor is ejected with a relative velocity opposing the mothership's orbital velocity, so the nominal release conditions of the meteor are given by: where the negative sign v rel indicates it is shot along the -X axis. Thus, the initial orbital velocity of the meteor with respect to an inertial, geocentric reference frame is easily obtained as the sum of the mothership's orbital velocity vector and rel , where the latter is expressed in coordinates of the already defined orbital frame, and as a function of the yaw angle (i.e. the rotation around the nadir) and pitch angle (i.e. the rotation around the second axis), as follows: Furthermore, in order to account for the mass variation during the reentry, this must be included within the state vector. The initial mass of the meteor depends on the initial radius R m0 and the initial density m0 through the following equation: A spherical shape is assumed with a uniform density distribution, such that R m0 = 5 mm and ρ m0 = 5000kgm −3 .

Formulation
The orbit propagation for the meteor problem has been carried out with Cowell's method [18]. It performs the numerical integration of the perturbed equations by using Cartesian coordinates, usually referred to an inertial frame, and characterized by six ordinary differential equations. The system is properly extended to describe the meteor motion by including an additional equation related with the mass loss: where is the Earth's gravitational parameter, v atm is the velocity relative to a corotating atmosphere defined in Eqs. (1.a) and (1.b), p is the perturbing acceleration due to the atmospheric drag introduced in Eq. (1.a), but in dimensionless form, namely: and [ , ] is the Cartesian state vector with respect to an inertial frame, expressed in dimensionless form using the following characteristic magnitudes to define dimensionless variables:

Design of the Monte Carlo Analysis
The calculations have been carried out in a computer with a 32 core Intel® Xeon® Gold 6130 processor @ 2.10 GHz under a 64-bit Ubuntu 20.04 LTS operating system. The solution provided by Cowell is integrated with Matlab's built-in ode113 integrator using a tolerance of 10 -13 . This value is fixed for both relative and absolute tolerance because all the equations in Cowell are dimensionless, and therefore it makes sense to fix the same value to both tolerances. Besides, in order to assess the influence of several variables in the propagation, and therefore in the expected requirements of the mission, a statistical analysis is performed by means of Monte Carlo (MC) analysis. The study is first performed by leaving aside the mutual interferences or cross-coupling that may exist between the different variables, i.e. varying only one variable at a time, while leaving the others frozen; this concession enables getting a clearer insight of the individual contributions of each variable of the sensitivity analysis. Later on, coupling effects for the release parameters are accounted for by means of two-and three-dimensional sensitivity analyses.

Nominal Mission
From our numerical simulations, the nominal values chosen for the experiment ensure the full mass depletion at an altitude of 70.7 km over the Earth ellipsoid (assuming the WGS84 model) and at a geodetic location of ( −146.7176 • West, 0.8558 • North) after a flight time of 15 min and 25 s.

Input Variables of the Experiment
The different variables involved in the definition of the mission will be classified in three different groups, since the variables within each of these groups needs to be modeled differently in the MC runs. These three groups are, respectively: 1. Engineering Design Parameters are those we can simply refer to as "design variables", in the sense that their value is a direct outcome of the experiment design, but once their value set, it is not subject to uncertainty, since their value can be measured and quantified prior to the experiment and remains constant thereafter, and thus these variables have known values. Within this category we consider the initial radius (R m0 ) and the initial density (ρ m0 ) of the meteor. Indeed, once their nominal values are defined, they can be measured on ground and during the experiment they are not subject to any uncertainty that requires a statistical modeling. Therefore, from a MC analysis viewpoint, it is most interesting to vary their value within a grid of prescribed, evenly spaced values for the sake of performing a parametric study within the domain of definition of these variables and provide the experiment designer with valuable insight of the consequences of varying these values.

Environmental and orbit-related variables
We sort within this category variables whose nominal or expected value can be estimated a priori, but is unknown or subject to uncertainty during the actual experiment; however, the uncertainty associated to these values is not necessarily relevant to be treated in a statistical way, either because their uncertainty is hard to model, they can exhibit unexpected an behavior, or because a statistical modeling is in principle not of interest. On the one side, we shall consider variables related to the space environment modeling, specifically related to the solar activity, such as the F 10.7 factor and the geomagnetic index Ap, because these variables exhibit a relative large uncertainty which is hard to match to a probability density function, so their values during the experiment can only be predicted roughly (although they can be measured in nearly-real time or in post-flight analysis); on the other side, we shall also include orbital variables such as the semi-major axis, because even though a nominal value is set on the design of the experiment, during the experiment, the actual value of this orbital parameter will exhibit some deviations due to the orbital perturbations. In the case of the aforementioned three variables, from a MC analysis viewpoint it is not interesting to approach their variations from a statistical modeling viewpoint, but instead using a grid of prescribed, evenly spaced values, for the sake of performing a sensitivity analysis and understanding how variations on these parameters may affect the outcome of the experiment. 3. Release conditions Within this category we shall include all variables related to the ejection phase of the meteor, namely the magnitude of the relative release velocity, v rel , and its direction, given by the pitch and yaw angles ( , ). These variables have associated uncertainties, since their values depend on either mechanical devices that cannot exactly reproduce the same ejection velocity, or the attitude of the spacecraft, which is also subject to a certain degree of uncertainty stemming from the attitude determination process. However, contrary to the previous two categories, it is now reasonable to model these variables statistically, assuming these are random variables and have an appropriate probability density function. To this end, normal distributions will be assumed for each of these three variables, centered at their nominal value and with a standard deviation that is consistent with assuming that 99% or the values are within prescribed intervals. As a result of treating these variables statistically, it becomes now possible to provide a statistical treatment also to the output variables of the experiment, as we shall see in the following section.

Output Variables of the Experiment
In order to assess the performance of the experiment, some "observables" or "output variables" of the experiment need to be monitored for each of the numerical simulations to be carried out in the following section. These will be the following: Time series for the state vector components and derived quantities, in particular the instantaneous mass and the visual magnitude of the meteor will be monitored both, as a function of time and as a function of the altitude. Also, the geodetic longitude, latitude and altitude will be monitored as a function of time.
Final values of the mass and the geodetic longitude, latitude and altitude will be recorded at the end of the experiment, either at the instant when the meteor fully depletes, or when it reaches (if it does) the Earth's surface.
If the meteor reaches the Earth's surface, the mass fraction that survives the reentry will be recorded. The minimum value of the time series for the visual magnitude during the reentry will be tracked too, which represents the maximum luminous intensity that the meteor reaches during the experiment. When Monte Carlo analyses are performed, dispersion plots and histograms will be used to display the distributions for some of the aforementioned observables.

Grid Sampling and Monte Carlo Analyses
Grid Sampling (GS) and Monte Carlo (MC) analyses will be performed in increasing order of complexity. First, one-dimensional analyses will be performed, where only a single variable will be varied while the rest are set to their nominal values, so the effects that varying each specific variable has upon the experiment results can be gauged and quantified. Afterwards, several two-dimensional Grid Sampling and Monte Carlo runs will be performed, where two are varied simultaneously, so their combined effect upon the experiment results and the output variables can be assessed. For these multi-dimensional analyses, two approaches have been employed: when combining two variables for which we don't have uncertainty information, and thus an evenly spaced grid is to be considered, a multi-dimensional GS is performed by constructing a 2D grid of sample points that stem as combinations of the prescribed grid value for each of the variables; when combining two variables for which we do have a statistical uncertainty model, we first combined their probability density functions to construct a multivariate probability density function using a Gaussian mixture model, so we can then sample as many combinations of input variables as we need and perform a MC run. Thus, in all cases, the procedure to set up a GS or MC run consists on defining: (a) the reference solution; (b) the variable(s) to be varied and the range where their values are to be varied; and (c) the number of grid points or random samples to be taken, and whether these are to be taken from a predefined grid of evenly spaced values, or from a univariate or multivariate random variable distribution, respectively.

Numerical Results
In this section the results of the different Grid Sampling (GS) and Monte Carlo (MC) runs are presented and explained. In the following figures, the blue line represents the reference or nominal solution, the gray lines represent each of the cases that the run has simulated either following a normal distribution or a grid of evenly spaced values, and the dashed lines (only visible in one-dimensional runs) represent the cases associated to the boundary values of the interval where the considered variable is varied. Histograms are used to visualize the distributions of the most relevant output variables by using either the end value, maximum or minimum value within the time series of each variable; in some cases, we will display histograms along with the best-fit Normal probability density function estimated from the output variables, so the validity of Gaussian hypothesis of the results can be tested. The final purpose of these GS and MC analyses is to assess whether the mission requirements are satisfied [19], to understand the dependence relations between input and output variables as well as their combined effects, and to assess the effects of a release different to the nominal value in any given input parameter. In order to keep an ordered approach, we shall investigate the effect of the input variables following the three groups or categories previously defined in Sect. 5.2.

Effect of the Engineering Design Parameters
In this section we study the effects of varying two of the input variables, namely the initial radius of the meteor, R m0 , and its density, m0 . Figure 4 shows how the visual magnitude and the mass of the meteor change along their trajectories as a function of the altitude over the Earth ellipsoid. To this end, the one-dimensional GS run accounts for a varying initial radius of the meteor ranging from 1mm to 10mm , where an evenly spaced distribution is assumed for the values of this input variable. In this figure, the light blue line in the solution profiles refers to the solutions stemming from the nominal conditions and parameter values; the gray lines represents each of the simulation cases considered in the sensitivity analysis; the dashed, black line indicate the envelope, determined by the parameter values farthest away from the nominal value; and the red dot indicates the nominal, initial value for the considered parameter. One of the immediate consequences of varying the initial 1 3 radius is that the initial meteor mass changes accordingly, so not only the altitude for complete depletion will vary, but there will also be cases where the meteor can actually reach the surface of the Earth before complete ablation. Therefore, it is also interesting to study how deep into the atmosphere the meteor can survive, and what mass fraction of the meteor makes it to the ground, as a function of the initial radius of the meteor; this is also illustrated in Fig. 4, where it can be observed that an initial radius of 6 mm ensures the meteor will reach the surface of the Earth, and beyond that threshold radius, an increasing amount of material will survive the reentry.
Another interesting aspect to observe is the brightness of the meteor during the reentry. The brightness, measured in mag0 units, is also displayed in Fig. 4, where a smaller value of the visual magnitude indicates a higher luminous intensity in the visual wavelength. Thus, it can be observed that the atmospheric entry conditions are such that the brightness peaks at about 80 to 120 km, depending on the initial radius. In particular, a higher meteor temperature (i.e. a higher velocity and higher atmospheric density) and a larger surface area will yield a more intense luminous emission. Consequently, the larger the meteor is, the deeper into the atmosphere that its luminous peak will occur, and the further away it will traverse along its orbit before either depletion or ground impact. In this regard, another interesting proxy to analyze is the geodetic distance that the meteor's final position (either at depletion altitude or on ground) will cover when projected onto the surface of the Earth's ellipsoid (i.e. the traversed groundtrack), depending on the initial size of the meteor. For the particular values considered in this MC run, the final depletion/impact covers a groundtrack distance spanning up to 442 km in the along-track direction, from (−146.4 Ultimately, the initial radius of the meteor affects its initial mass, and so does the meteor's density, m0 , which can be varied by choosing an alternative material for the meteor. Thus, varying the meteor's density has qualitatively the same effects on the experiments output variables, and therefore we opt to omit any additional figures for the sake of concision. For the density we chose to vary its values within the range of 1,000 to 10,000 km m −3 using an evenly spaced distribution, which would yield a final depletion/impact groundtrack distance covering up to 277 km in the along-track direction, from (−146.6 • West, 3.497 • North) to (−146.8 • West, −0.8073 • South).
However, increasing the density augments the actual mass of the meteor while maintaining its surface area constant, as opposed to varying the initial radius, which yields changes in both, the mass and size of the meteor. Therefore, by conveniently varying both, R m0 and m0 , it should be possible for the experiment designer to achieve any desired combination of initial mass and surface area of the meteor so that either the depletion altitude or the mass fraction at ground impact can be chosen at will; alternatively, it is also possible to design the maximum luminous intensity that the meteor reaches, as well as the geodetic altitude at which it occurs. Therefore, it is interesting to analyze the combined effects of varying these two design parameters simultaneously by means of a two-dimensional GS where both variables are varied so they take combinations of the aforementioned values. This is illustrated in Fig. 5, where the variations of the visual magnitude and actual mass of the meteor as a function of time exhibit qualitatively the same behavior as in Fig. 4, only that now the dispersion of the results is larger, because the considered combinations of the input values cover a wider range of initial masses, ranging from the smallest and most lightweight meteor of 0.042 g (radius of 1 mm and density of 1,000 kg m −3 ) and up to the largest

Effect of the Environmental and Orbit-Related Variables
In this section we shall study the consequences of varying three different input variables: on the one side, two variables related to the solar activity, namely the F 10.7 solar flux and the geomagnetic index Ap, which directly affect the atmospheric density; on the other side, one variable related to the mothership's orbit (and therefore the release position of the meteor), in particular its semi-major axis, a.
The solar activity is of paramount importance for any orbit propagation that accounts for the drag perturbation, and in particular an orbital reentry. The two The Journal of the Astronautical Sciences (2022) 69:1347-1374 key parameters that model the solar activity are the F 10.7 solar flux and the geomagnetic index Ap. We shall first look into the former by means of a one-dimensional GS assuming a range of evenly spaced values for the solar flux ranging from 60 to 240, so they correlate with a minimum and maximum solar activity, respectively. Figure 6 displays the visual magnitude and the mass of the meteor change during the reentry trajectory as a function of the altitude. As opposed to the results of Sect. 6.1, varying the solar activity has in practice a negligible effect upon the depletion altitude or its geodetic position; indeed, in all the considered cases the meteor fully depletes within an altitude range of barely 250 m regardless of the solar activity, and the maximum groundtrack covered by the final depletion point spans only across 1.69 km in the along-track direction, from (−146.7 • West, 0.8675 • North) to (−146.7 • West, 0.841 • North) . Thus, the only meaningful effect of the solar activity seems to be in the rate of change of the mass during reentry. A higher solar activity rises the atmospheric density and temperature at the upper layers of the atmosphere, which in turn raises the ablation rate at higher altitudes and raises the luminous intensity of the meteor through the upper atmosphere; however, neither the full-depletion altitude/position, or the altitude and intensity where the visual magnitude peak occurs, exhibit any meaningful variations and they remain in practice quite invariant.
The geomagnetic index was also varied in a one-dimensional GS with evenly spaced values ranging from 1 to 40, which correlate with a low and high geomagnetic activity, respectively. Their impact upon the observed output variables was similar to that observed in Fig. 6, only that the magnitude of the variations were far less pronounced to the point that they were barely noticeable in the plots, so we chose to omit these figures. Hence, the combined effect of simultaneously varying the F 10.7 solar flux and geomagnetic index Ap is dominated by the value of the solar flux, and thus results look like those of Fig. 6. Consequently, we arrive at the conclusion that the solar activity does not affect the experiment results in any meaningful way, as the altitude of full depletion exhibits only marginal variations, as shown in Fig. 7. The other variable to be analyzed in this section is the mothership's orbit at the time of releasing the meteor, in particular its semi-major axis, a. Even though a nominal orbit will be defined for the mothership, it is important to be aware that the orbital perturbations can impose short terms variations upon its value, so at the time of release the actual orbital altitude of the meteor may be slightly different from the intended value, and thus it is important to gauge how this may affect the experiment results. To this end, we shall perform a one-dimensional Grid Sampling assuming an evenly spaced set of semi-major axis values ranging 10 km above and below the nominal value.
The main effects of varying the semi-major axis of the mothership are twofold: on the one side, the meteor is released at a slightly lower altitude; consequently, since the mothership is at a lower orbit (thus a lower period) it is released with a slightly larger initial geocentric orbital velocity. The combination of these two aspects has a series of consequences that are illustrated in Fig. 8. Interestingly enough, the altitude of full-depletion does not vary significantly, yet meteors released at a lower altitude reach the Earth's atmosphere earlier and therefore have shorter flight times, whereas meteors released at higher altitudes have longer flight times. Interestingly, for the considered small variations, the release altitude and the flight time correlate nearly linearly, and consequently so does the subsatellite position where the full depletion occurs, as evidenced by the uniform spacing of the geodetic positions where the complete ablation occurs. Additionally, when the visual magnitude and actual mass are plotted as a function of the altitude, they all overlap to the point that, on the plot, they become visually indistinguishable from one another, so these magnitudes must be plotted as a function of time instead; it is like this that one can observe that meteors released at a higher altitude have longer flight times mainly because they reach the upper layers of the atmosphere later, and thus the ablation is delayed compared to meteors released from a lower altitude. This effect also delays (for meteors released at higher altitude) the moment where the minimum visual magnitude (maximum luminous intensity) is achieved, although it does not change its peak value.

Effect of the Release Conditions
In this section the consequences of varying the meteor's release velocity will be analyzed. In particular, the impact of the following three input variables will be analyzed: the magnitude of the relative release velocity, v rel , and the yaw ( ) and pitch ( ) angles, which provide the direction in which the meteor is ejected from the mothership. These three variables have associated uncertainties which are most conveniently modeled under a statistical approach; thus, these variables will be treated as random variables with associated probability density functions modeled as normal distributions centered at the nominal values, and with a standard deviation that is compliant with a desired interval width for the values that these variables can take with a given confidence level. Thus, we shall define the following probability density functions: The values of the standard deviations are chosen so that, for a large number of samples, 99% of the v rel samples fall within its nominal value ±3% , and so that for the yaw and pitch angles, 99% of the samples fall within an interval of ±3 • , which is thought to be consistent with the expected accuracy with which the mothership can determine its own attitude state. These values are compatible with the capabilities of the release system, as described in [20]. Figure 9 represents the results of a one-dimensional MC run for the v rel variable based on 1,000 random samples. The main effects of varying the relative release velocity, and thus the initial orbital velocity of the of the meteor, are similar to those observed in Fig. 8 for varying the mothership's semi-major axis, in the sense that a lower relative release speed (hence, a larger initial orbital velocity) slightly increases the flight time by delaying the atmospheric reentry, and thus allowing the meteor to cover more groundtrack before its full depletion. The main difference resides in that here the release occurs always at the same orbital position, and only the orbital velocity is varied. Note that for the considered uncertainty model, variations in v rel = N( = 0 • , 2 = 1.16279 deg 2 ) Interestingly enough, small variations of the release velocity also yield a nearly linear relation with the flight time, and therefore the full-depletion position in the along-track direction. This means that the flight time and the geodetic longitude, latitude and altitude are all mutually correlated with the initial release velocity; consequently, the assumed normal distribution for v rel should translate into normal distributions for the aforementioned output variables too. This is illustrated in Fig. 10, where best fit normal distributions have been overlaid on top of the histograms showing the occurrences of the geodetic longitude, latitude and altitude. It can be observed that for all these variables the normal distributions exhibit a good fit with the observed values. In particular, the final geodetic longitude at the instant of full depletion can be modeled by a normal distribution N( = −146.7171 • , 2 = 0.022699 deg 2 ) , for the latitude we get N( = 0.86472 • , 2 = 0.39207 deg 2 ) and for the altitude N( = 71.6997 km, 2 = 0.20162 km 2 ) . Note that the fit of the altitude distribution in Fig. 10c can be improved assuming a bi-modal distribution instead, to better capture the higher altitude cases. A one-dimensional MC run with 1,000 samples varying the pitch angle, , yields qualitatively the same behavior as in Fig. 9, as varying the in-plane orientation of the relative release velocity vector, rel , is equally equivalent to a change in the initial geocentric orbital velocity of the meteor, and thus we shall not include additional figures. For all the considered cases the meteor still fully depleted within a range of geocentric altitudes spanning 464 m while the flight times diverted from one another by a maximum of 62 s and the peaks of the visual magnitude remained of the same intensity. The maximum groundtrack covered by the final depletion point spans across 252 km in the along-track direction, from (−146.6 • West, 2.797 • North) to (−146.8 • West, −1.108 • South) . Attending at the histograms of Fig. 11 we observe that the distributions of the geodetic longitude and latitude still follow a normal distribution, in this case N( = −146.7194 • , 2 = 0.03755 deg 2 ) and N( = 0.82784 • , 2 = 0.57397 deg 2 ) respectively, but the altitude now no longer resembles a normal distribution, but it seems to obey instead a bi-modal distribution.
A one-dimensional MC run with 1,000 samples varying the yaw angle, , however, yields a different structure of the output variables compared to the previous two cases. The effect of the yaw angle if that of providing the initial orbital velocity of the meteor with an off-plane component, which directly affects both, the geodesic longitude and latitude of the meteor's full depletion point providing it a symmetrical structure with respect to the orbital plane, as can be observed in Fig. 12. Note that the net effect of this off-plane velocity component is comparatively much smaller than in-plane variations in terms of modifying the orbital geometry; consequently, note that the variations upon output variables such as the geodesic longitude, latitude, altitude, mass variation and visual magnitude, are so subtle that they are barely noticeable in a plot unless we zoom in, so in practice they can be considered independent from the yaw angle for the considered range of values. In fact, in all cases the meteors fully depleted within a range of altitudes spanning 45 m, the flight times diverted from one another by no more than 1.5 s and the peaks of the visual magnitude remained of the same intensity. On the contrary, the final geodetic position at full depletion can exhibit a comparatively wider dispersion on the ground map, of up to 36 km in longitude ( 0.3235 • ) and 10.46 km in latitude ( 0.09396 • ). Attending at the histograms, only the longitude at final depletion exhibits a normal distribution (with parameters For the input variables that determine the release conditions of the meteor, it is extremely interesting to perform two-dimensional MC runs to study the combined effect that these variables have upon the observables or output variables of the experiment. Since the release velocity and pitch angle have similar effects and can only produce a dispersion of the final depletion time in the along-track direction, it is most interesting to analyze their combinations with the yaw angle, since the latter is the only of these three variable capable of providing a noticeable dispersion in the across-track direction. To this end, a bivariate normal distribution was constructed based on pairs of the same and 2 values previously used when individually analyzing these variables, from which 5,000 samples were used for the MC run. Results illustrated in Fig. 13 for the combined effect of simultaneously varying the pitch and yaw angles. The left plot represents the dispersion of the on-ground projection of the final positions of the meteor at time of full-depletion or ground impact, along with the "minimum volume enclosing ellipse", which represents the "ground error ellipse", in dotted, black line. The remaining ellipses represent contour lines of the best-fit Gaussian mixture distribution model representing the occurrences for the final longitude and latitude at full-depletion or ground impact; these levels correspond to non-normalized values [1,2,3,4,5,6], respectively, 1 corresponding to the outermost contour level, and 6 to the innermost. The red circle indicates the depletion point of the nominal meteor trajectory, while the black asterisk pinpoints the geometric center of the enclosing ellipse. For the combined variations of the pitch and yaw angles, all cases still lead to a full atmospheric ablation within altitudes spanning 493 m (ranging from 71.60 to 72.083 km) and flight times spanning from 14.94 to 15.98 min (i.e., no more than 62.5 s in difference), and the error ellipse has semi-major axis of 238.086 km ( 2.1388 • ) and a semi-minor axis of 2.1388 km ( 0.16818 • ). Thus, the distributions of the geodetic longitude and latitude follow normal distributions N( = −146.7185 • , 2 = 0.055579 deg 2 ) and N( = 0.85258 • , 2 = 0.52952 deg 2 ) respectively, whereas the altitude exhibits a bi-modal distribution as in Fig. 11, as they are dominated by variations in the pitch angle. For the combined variations of v rel and the outcome would have been similar, both qualitatively and quantitatively, thus showing that the experiment is robust against the expected uncertainties in the release conditions, so we decided not to include additional figures for the sake of concision.
Finally, we have considered an interesting addition to include a hypothetical situation where the experiment goes widely off the nominal values. To create such a scenario, we intentionally increased the uncertainty of the variable v rel from ±3% of its nominal value to ±35% , in order to assess the potential consequences. In such a hypothetical situation, the dispersion of the final meteor positions would increase substantially (thus enlarging the enclosing ellipse on the ground map), to the point where some of the meteor samples would even survive the reentry and reach the ground, as shown in Fig. 14. However, even then, the probability of the meteor reaching the Earth surface would be as low as 1.2% based on this MC run, and from all the samples considered the largest mass fraction reaching the Earth was of

Sensitivity Analysis
Finally, the quantitative effects of variations upon the nominal values of the mission parameters, and how these produce changes in the observables of the experiment, can be highlighted by means of a classical sensitivity analysis plot, as the one shown in Fig. 15, which showcases the relationship between variations in the input variables of the experiment, and the resulting variations in the output variables. The subplots of Fig. 15 illustrate how prescribed variations on the input variables yield changes in four observables of the experiment, namely: the geographic longitude, latitude and altitude of the meteor at the time of full-depletion, and the peak value of the visual magnitude during the experiment. The purpose of these sensitivity plots is to allow for the relative quantification of the experiment outcomes under comparable deviations in each of the input variables or parameters, thus to highlight the relative importance of each input variable towards the success of the experiment, and therefore to point out how sensitive the experiment is, in quantitative terms, to potential deviations from the nominal values of these input variables. To this end, variations of ±1% have been considered for the input variables R m0 , m0 and v rel , as well as for the initial orbital altitude of 375 km ; for the F 10.7 solar flux, more realistic variations of 10% were considered; for the geomagnetic index, A p , variations of ±1 were considered, since by definition this index can only take integer values; for  Fig. 15, whereas negative variations are indicated by red bars, which allows to visualize also whether an increase/decrease in an input variable correlates with an increase or a decrease in a certain output variables. As already concluded in previous sections, these results confirm that the experiment results are significantly more sensitive to the values of the engineering design variables and the release conditions, whereas orbital altitude (alternatively, the semi-major axis) has a comparatively lower impact, and the environmental variables, namely the F10.7 solar flux and geomagnetic index, have a comparatively little impact in the experiment outcomes; the latter is fortunate, since the actual values of these environmental variables during the experiment are hard to estimate with precision beforehand. Regarding the design variables, R m0 and m0 , although they can have a significant impact on the experiment, deviation in these variables are also easy to minimize, since they can be achieved with great accuracy. Therefore, it is the release variables that are most sensitive in practice, since their actual values during the experiment are subject to operational errors, as well as to uncertainties in the attitude determination of the spacecraft.

Achievements
This paper presents the recent efforts to improve the characterization of the meteors' trajectory and safe ablation. The dynamics during the atmospheric entry are challenging to simulate because the model presents a large amount of uncertainties. Consequently, the meteor trajectory was assessed through a statistical analysis of the parameters involved in the physical modeling by considering Cowell's special perturbation method. The results confirm that the largest influence in the trajectory of the artificial meteor comes from the engineering design parameters and the ejection parameters. They can not only change the final position of the meteor but also the final state and consequently the corresponding altitude of full depletion. The engineering design parameters have the advantage that they can be chosen a priori with very good accuracy and, therefore, the problem can be avoided by selecting appropriate values for the mission. Special attention has to be given to the ejection parameters because they are the only ones which present uncertainties coming from different sources, like the attitude measurement and determination. However, the environmental parameters are not as relevant to the trajectory, which becomes a positive aspect taking into account that the exact value of these parameters cannot be determined a priori. Accordingly, we can conclude from our analysis that small deviations from the nominal values are still compliant with the requirements of the mission [19].

Future Work
In this paper the physical model has only considered the mass ablation. Ongoing academic works include, but are not limited to: obtain more realistic physical and dynamical equations by exploring the effects of additional mechanisms like the melting and evaporation of the material [21],modeling of ablation within metals, ceramics, meteors and ALE materials; coupling of material ablation with flow aerodynamics; computation of the spectral properties of ALE materials. In a future work, the analysis will also be refined by including a comprehensive modeling of the particle shape deformation and its effect on the aerodynamic coefficients and trajectory.