Age of Geminids Derived from the Statistics of Meteoroid Orbits

Statistical analysis of samples of the orbits of celestial bodies is complicated by the fact that the Keplerian orbit is a multidimensional object, the coordinate representation of which nonlinearly depends on the choice of orbital elements. In this work, using the construction of the Fr\'echet mean, concepts of mean orbit and dispersion of the orbit family are introduced, consistent with the distance function on the orbit set. The introduced statistical characteristics serve as analogs of sample mean and variance of a one-dimensional random variable. Exact formulas for calculating the elements of mean orbits and dispersion quantities with respect to two metrics on the orbit space are derived. For a large sample of meteoroid orbits from the Geminid stream, numerical simulations of orbit evolution over 20,000 years in the past are conducted. By analyzing the dependency of statistical characteristics on time, estimates for the age of the stream and the gas outflow velocity are obtained under the assumption of the birth of the Geminids due to the rapid destruction of the cometary nucleus.


Introduction
In the works Kholshevnikov (2008) and Kholshevnikov et al (2016), a distance function on the space of Keplerian orbits was constructed, directly reflecting the difference in the values of the integrals of the two-body problem on solutions.Specifically, if x and y are two non-rectilinear orbits, the distance between them is defined as where u x , v x and u y , v y are pairs of orthogonal vectors aligned with the angular momentum vector h and eccentricity vector e corresponding to orbits x and y, given by the formulas where κ 2 is the gravitational parameter, and r and ṙ are the position vector and velocity of the point.The metric on the orbit set is a valuable tool for studying the genealogy of small bodies in the Solar System, as it allows characterizing the proximity of two points on the five-dimensional manifold of osculating Keplerian orbits, which is nontrivially topologically structured.
In the context of meteoroid streams, the metric (1) is used both to search for parent bodies (Kokhirova et al (2018); Sergienko et al (2020)) and to determine the stream to which a specific body belongs (Andrade et al (2023)).
In this work, we will define and provide methods for calculating statistical characteristics of the orbit sample: mean and dispersion, analogous in meaning and properties to the sample mean and variance of a random variable.To achieve this, we will apply the construction of the Fréchet mean in the metric space (Fréchet (1948)) to the orbit space with the metric ϱ 2 .
Using the dispersion, reflecting the "size" of the stream in terms of its concentration around the mean value, we will estimate the age of the Geminid meteoroid stream and the gas ejection velocity, assuming its formation as a result of the rapid destruction of the cometary nucleus.
Our simulation of the evolution of the orbits of meteoroids of the stream over 20 000 years into the past shows that the dispersion of the sample increases over large time intervals, due to perturbations.The ensemble of meteoroid orbits constituting the stream "diffuses" over time, merging with the background.In accordance with this, we propose to search for the age of the Geminids among the minima of the dispersion of the orbit sample close to the present time.The success of this approach is not guaranteed in general since any orbit sample obtained from meteor observations represents only a small part of the variety of meteoroid stream orbits.However, as will be seen later, in the case of the Geminids, a reasonable estimate of the stream's age can be achieved.

Mean in the space of Keplerian orbits
The operation of calculating the mean for a sample of Keplerian orbits does not have a natural and unambiguous definition.Often, for this purpose, the procedure of calculating the arithmetic mean of orbit elements is used.However, such a method has a drawback: in different systems of elements, the mean orbits of the same family of bodies may, in general, differ.The reason for this lies in the nonlinearity of the transformations between element systems.
This drawback of element-wise averaging of orbits was noted in the article Jopek et al (2006).It was suggested to calculate the mean orbit as the best least squares approximation of orbits in the family, parameterized by the vector (h, e, E), where h and e are the angular momentum and eccentricity vectors, and E is the total energy of a unit-mass body on the orbit.The article provides a numerical algorithm for solving the minimization problem, using the arithmetic means h, e, and E as initial approximations.The question of the uniqueness of the solution was not considered.
This approach can be generalized as follows.Starting from the work Southworth and Hawkins (1963), any study of a family of bodies with close orbits directly or indirectly uses some numerical criterion of orbit proximity.Let it be denoted by D. As long as, in the researcher's opinion, D adequately reflects the concept of proximity, it is natural to define the mean based on this criterion.Specifically, for a family of orbits x k , k = 1, . . ., n, we will call the mean orbit x the one that provides the global minimum of the function The set of points of global minimum of S is known as the Fréchet mean (Fréchet (1948)), when D is a metric on the set of all orbits.This type of mean represents a deeply meaningful generalization of the arithmetic mean to an abstract metric space.For the Fréchet mean in an arbitrary space, an analogue of the strong law of large numbers holds (Ziezold (1977); Molchanov and Molchanov (2005)), and in the case of a Riemannian manifold, uniqueness conditions are derived (Bhattacharya et al (2003)), along with a generalization of the central limit theorem (Bhattacharya and Patrangenaru (2005)).
The element-wise mean orbit mentioned above falls within the definition of the Fréchet mean if we set where ϶ 1 , . . ., ϶ 5 are the elements of the orbit in the chosen system.The mean from Jopek et al (2006) corresponds to the metric The value of the function (3) at the Fréchet mean, that is, the minimum over the entire orbit space H, represents a natural generalization of the sample variance, describing the characteristic "size" of a random sample in terms of the root mean square deviation from the mean.Calculating the Fréchet mean, that is, minimizing the function (3), is usually a challenging task, but for the metric ϱ 2 , it is possible to obtain exact formulas.

Mean with respect to the ϱ 2 metric
From definitions (3) and (1), it follows that to find the Fréchet mean in the ϱ 2 metric for a family of orbits defined by vectors (u k , v k ), k = 1, . . ., n, one needs to minimize the quadratic form defined in R 6 .Vectors u and v define a curvilinear orbit only if u ̸ = 0 and Thus, the form (5) should be minimized subject to the condition (6).The solution to this problem, as shown in Appendix 5, yields the following values for the vectors u and v of the mean orbit: where In exceptional cases where ū is collinear with v, µ 2 calculated by formulas (8) equals 1.In these cases (more details are discussed in Appendix 5), either there is no solution or there are infinitely many solutions.However, the collinearity of averaged angular momentum and eccentricity vectors indicates an extremely high degree of dispersion in the sample.
In problems of averaging samples of meteoroid orbits, such a situation is unlikely.For a family of orbits that are close to each other, the magnitude of |µ| is small.Among meteoroid streams from the Jenniskens et al (2016) catalogue, its maximum value is 0.028 (stream 164 Northern June Aquilids), and the median is 0.0005.A first-order approximation in terms of µ, easily obtained from (7), provides an insight into the deviation of the parameters u, v of the mean orbit from the arithmetic means of the same parameters of the family.The value of the function S on the mean orbit naturally generalizes the concept of the root mean square deviation, characterizing the scatter of orbits in the family.Let's denote this value as S 2 and express it in terms of the elements of the family: where The relationship between the Keplerian orbit elements and the components of vectors u and v is determined by equalities (Kholshevnikov et al (2016)) where e, p, i, ω, Ω are the eccentricity, semi-latus rectum, inclination, argument of pericenter, and longitude of the ascending node of the orbit, respectively.It is not difficult to obtain formulas for the inverse transition, from u and v to e, p, i, ω, Ω.

Mean with respect to to the ϱ 5 metric
When studying families of orbits that have undergone long-term evolution, the metric ϱ 2 may be too subtle a tool.Under the influence of small perturbations, the argument of pericenter and longitude of the ascending node of an orbit can change rapidly compared to the other elements.A metric that neglects the difference in rapidly changing orbit parameters is defined in Kholshevnikov et al (2016) by the formula The distance ϱ 5 between orbits x and y is equal to the minimum distance ϱ 2 between orbits X and Y, all Keplerian elements of which, except for Ω and ω, coincide with the corresponding elements of x and y.
Obviously, ϱ 5 (x, y) = 0 if x and y differ only in the argument of pericenter and longitude of the ascending node.Thus, the orbit space H is divided into equivalence classes: for any two elements of the same class, ϱ 5 equals zero, and formula (9) defines the distance between the classes.In Kholshevnikov et al (2016), it is shown that ϱ 5 satisfies the triangle inequality.That is, the described set of equivalence classes, equipped with the distance ϱ 5 , forms a metric space.
From the results of Milanov (2018), it follows that this space isometrically embeds into Euclidean space R 3 .The image under this mapping is a convex set with a cut-off ray, determined by the requirement u ̸ = 0.This embedding allows reducing the problem of computing the mean orbit to finding the ordinary arithmetic mean of vectors in Euclidean space.The following formulas for computing Keplerian elements of the mean orbit were obtained in Milanov (2019).We present them below.
where p j , e j , i j , j = 1, . . ., n are the semi-latus recta, eccentricities, and inclinations of the orbits of the averaging family of size n, and p, e, i are the corresponding elements of the sought mean orbit.Since the mean in the considered space is a unique image of the arithmetic mean in Euclidean space, the obtained solutions are unique.The mean defined by formulas (10) exists if the plane of at least one orbit does not coincide with the reference plane or u k ̸ = 0.This limitation is a consequence of the requirement u ̸ = 0 for the mean orbit.
The dispersion, S 2 5 , defined by the general formula (4), is expressed in Keplerian elements as follows: 3 The dispersion and the age of the Geminid meteor stream The difficulty of statistical analysis of a sample of Keplerian orbits lies primarily in the fact that an orbit (without considering the body's position on it) is a point in a five-dimensional manifold.The introduced quantities S 2 and S 5 simplify the task, providing numerical characteristics of the orbit sample analogous to the dispersion of a one-dimensional random variable.
It can be verified, by repeating the classical reasoning, that for the Fréchet mean µ and the quantity S defined by formula (4), Chebyshev's inequality holds in the form where P is the counting measure on the sample.Thus, as in the one-dimensional case, the purpose of the dispersion is to describe the characteristic size of the sample in terms of the probability of deviation from the mean.Next, we assume that for the ensemble of meteoroid orbits, this size was minimal at the time of the stream's formation.In other words, the age of the stream should be sought among the minima of the dependencies of S 2 and S 5 on time.We built these dependencies for the Geminid stream, modeling the evolution of the orbits of several thousand meteoroids over the past 20 000 years.Certainly, no matter how large a sample we take, we remain limited to orbits obtained from meteor observations on Earth's surface, that is, we "see" only a small part of the stream.To confirm the value of the age determined as the minimum of S 5 , we calculated the moments of the minimal distance between the orbit of the Phaethon which is the presumed parent body of the stream and the mean orbit of our sample.Finally, assuming a cometary scenario for the formation of the Geminids and limiting the gas escape velocity, we obtained an upper estimate of the age of the stream.

Geminid meteor stream
The Geminids are one of the most intense meteor showers, with the number of observed meteors per hour exceeding 100 at the peak of activity (Rendtel (2014)).Observations of the meteor stream have been documented at least since the second half of the 19th century (Greg (1872), Sawyer (1880)), and are currently conducted both from the Earth's surface (photographic, radar, video observations) and in space conditions, such as on the Parker Solar Probe (Cukier and Szalay (2023)).A notable feature of this stream, in addition to its significant intensity, is the nature of the parent object: in several studies, it is considered that the parent body is an asteroid, (3200) Phaethon, rather than a comet (Neslušan (2015), Yu et al (2019), Cukier and Szalay (2023)).
The Geminids are part of a larger complex of Phaethon-Geminid Stream Complex (PGC) together with the Daytime Sextantids, the study of which using ground-based radar and optics shows that the associated parent asteroid (2005) UD may have originated from the catastrophic destruction of a larger object that also gave rise to (3200) Phaethon (Kipreos et al (2022)).Nevertheless, in some studies, the connection between (2005) UD and the complex is disputed (Ryabova et al (2019)).Studies of the reflected spectrum from (3200) Phaethon and ( 2005) UD emissions show significant differences in the surface composition of the objects, which may indicate either the absence of any connection between the objects or the different effects of space weathering on the objects (Kareta et al (2021)).
Cukier and Szalay (2023) modeled the Geminid meteor stream and compared the results with data on the distribution of dust particles obtained by Parker Solar Probe; the core of the stream was found to be located outside the orbit of the asteroid (3200) Phaethon, which is more consistent with the stream originating from catastrophic disruption, rather than comet-like activity.Nevertheless, comet-like activity of the asteroid was observed during perihelion passages in 2009, 2012, and 2016, and there are models of the structure of the asteroid allowing for the formation of the complex by cometary activity.In this case, the thickness of the dust crust above the water-ice layer should not have exceeded 1 meter about 1000 years ago, and in the modern era, the activity is assumed in the vicinity of perihelion (Yu et al (2019)).A cometary scenario for the formation of the Geminids is also indicated by Ryabova et al (2019), without excluding the possibility of a scenario involving rapid ejection of volatile particles due to a catastrophic process.
Estimates of the age of the stream obtained by different researchers range from hundreds of years to tens of thousands.Based on the statistical difference in the orbits of the stream corresponding to meteoroids of different masses, Babadzhanov and Obrubov (1984) provides estimates in the range of 1600 − 19 000 years, depending on the particle density.Using numerical integration, Jones (1985) obtained a value of 6000 years.In Gustafson (1989), an interval of 1000 − 2000 years is given, and it is noted that with more significant influence of non-gravitational forces, the estimate will decrease.The value of 6000 years is mentioned in Rendtel (2005).By comparing the stream model with observational data, Ryabova (1999) finds the best agreement for an age of 2000 years.

Meteoroid orbits selection and the numerical integration
We used meteoroid orbit data obtained over several years of observations from the Global Meteor Network (GMN) (Vida et al (2021(Vida et al ( , 2020))), available on the GMN website (https: //globalmeteornetwork.org/data/traj_summary_data/).The selection of orbits for modeling the evolution was done yearly, covering the period from 2019 to 2023, resulting in five separate datasets.Particles identified with the Geminids meteor shower by GMN's internal algorithms were considered.In total, we investigated five datasets, comprising more than 53 thousands meteoroid orbits: 2116 for 2019, 5959 for 2020, 9968 for 2021, 15 795 for 2022, and 19 653 for 2023.
For each orbit, the meteoroid's position was recalculated to the Julian date ten days before the observation started in the respective year and the mean anomaly was adjusted accordingly.Then, the meteoroid's motion was computed backward in time for 20 000 years, with a time step of 0.005 years.
The integration of motion equations was performed using the numerical integrator collo (https://github.com/shvak/collo),developed by V.Sh.Shaidulin.Based on this integrator, a software package was developed to study the dynamics of various meteoroid streams represented in the Global Meteor Network.The package also allows for computing the average orbit, root mean square deviation, and conducting comparisons between the average orbit and the orbit of the assumed parent body in various metrics.
Text files with the obtained orbital elements for the entire simulation period in 5year increments are available at https://disk.yandex.ru/d/HLWHRZzcQwLVwQ. Exactly these data were utilized in the current study.
Orbits that experienced significant perturbations during the modeling period were excluded from further analysis.An orbit was considered strongly perturbed if the absolute value of the change in one of its elements during a five-year period exceeded the upper 2% quantile of all such changes for the respective dataset.Table 1 presents the threshold values for the absolute differences in elements, and Figure 1 shows histograms of the 99% absolute values of differences for the 2020 dataset, indicating the cutoff threshold.This filtering procedure excludes 5% to 6% of the orbits from each year's dataset.
The results described further were obtained for the 2020 dataset.We opted not Year ∆a (a.u.) ∆e ∆i

Mean orbit of the sample
The evolution of Keplerian elements of the mean orbit, calculated using formulas (7), is shown in Figure 2. A significant variation in inclination is noticeable, ranging from 13 • to 45 • , while the eccentricity changes within the range of 0.81 to 0.93.The semi-major axis undergoes changes of no more than 0.05 AU until approximately 15 000 years, after which this parameter sharply increases.The argument of pericenter, and longitude of the ascending node change more significantly than other elements, completing almost a full revolution over the simulated time span.The rapid changes in ω and Ω are accompanied by a larger spread of values compared to other elements, as reflected in the graphs depicting the standard deviations of orbital elements over time shown in Figure 3. Here, the sample variance is calculated for points on the unit circle with arguments i, arcsin e, Ω, and ω.The "angular" orbit element arcsin e is used to ensure that all compared quantities are angles.At an angle arcsin e, a segment connecting the center with the focus is visible from the vertex of the ellipse.The dispersion of Ω and ω significantly exceeds the variance of other elements.Therefore, along with the mean orbit with respect to the ϱ 2 metric, we will also consider the mean orbit in the ϱ 5 metric and the associated dispersion value S 5 defined by formula (11).
Recall that ϱ 5 is a pseudometric on the space of orbits that "glues together" orbits that differ only in the values of the longitude of the ascending node and the argument of the pericenter.The mean orbits with respect to the ϱ 2 and ϱ 5 metrics are close to each other for most of the simulation period.The distance ϱ 2 between them does not exceed 0.005 for the first 15 000 years and then sharply increases.

Dispersion of the orbits sample
Let's examine the evolution of quantities S 2 and S 5 over time.Defined by the general formula (4), these values represent analogs of sample variance on the manifold of orbits equipped with metrics ϱ 2 and ϱ 5 .Figures 4 depict the evolution of S 2 and S 5 over the entire simulation period and during the first 4000 years, which is of particular interest.Additionally, for the 2020 sample, orbit modeling was extended 10 000 years into the The intervals of monotonicity for the quantity S 2 closely resemble the corresponding pattern for S 5 .However, the minimum of S 2 is located later -around the 300-year mark.The position of this point is determined by a combination of decreasing dispersion for slowly changing elements p, e, i and rapid growth of the variance for Ω and ω (see Figure 3).We consider this minimum to be an excessively biased estimate of the age, with high uncertainty due to the low diversity of the sample in terms of Ω and ω compared to other elements.

Distance between the mean orbit and Phaethon's orbit
Alongside the meteoroid orbits, the orbit of Phaethon, the presumed parent body of the stream, was modeled over the same time span.Plots of distances ϱ 2 and ϱ 5 between the mean stream orbit and the asteroid's orbit are presented in Figure 5.The values of the distance functions, drawn in black, experience frequent small fluctuations, caused by oscillations of the elements of Phaeton's orbit.The colored curves depict the moving averages with a window of ±250 years.The global minima of the averaged dependencies are located at 1700 years for ϱ 2 and 1200 years for ϱ 5 .Thus, according to the modeling results, during the interval T 2 = [1200, 1700] years ago, the considered stream of orbits had the smallest size in terms of the root mean square deviation S 5 , and its mean orbit was closest to the orbit of the presumed parent body.
The distance ϱ 5 also exhibits a deep local minimum at 10 000 years, but the nearest minima of S 2 and S 5 differ by more than 2000 years.

Dispersion of the orbits sample and the ejection velocity
It is interesting to compare the obtained values of S 2 with the maximum possible values at the time of stream formation.Such a comparison will provide, on the one hand, an age constraint, ruling out obviously impossible scenarios, and on the other hand, a lower limit on the ejection velocity in a cometary scenario of stream origin.Assuming the Geminids are the result of a gas dust outburst from Phaethon occurring near the perihelion of the asteroid's orbit, we can estimate the dispersion from above with the maximum distance ϱ 2 between Phaethon's orbit and the ejected particle's orbit.
In terms of ejection parameters, this distance depends on the position of the parent body on its orbit and the relative velocity vector of the meteoroid.Estimates of the ejection velocity magnitude V calculated in Ryabova (2013) give values of 500 to 830 m/s for dust particles and 945 to 1368 m/s for gas, depending on the theoretical ejection model.The same work notes that these values are not sufficient to explain the observed width of the stream in ecliptic longitude, suggesting that the velocity could have been significantly higher.
Let t denote a moment in the past.Knowing the elements a, e, i, Ω, and ω of Phaethon's orbit at time t and additionally specifying the true anomaly θ, we can determine the point in space r(t, θ) where the asteroid was located at a time close to t with a given true anomaly and its velocity ṙ(t, θ).We introduce the quantity R 2 = R 2 (t, θ, V ), equal to the maximum distance between Phaethon's orbit and the orbit of a particle passing through the point r(t, θ), such that the particle's velocity at this point differs from ṙ(t, θ) by no more than V in magnitude.From the definition (4) of dispersion, it follows that at the moment of ejection, for any sample of meteoroid orbits X of size n: assuming that the particle velocities do not exceed V .The value of R 2 for given t, θ, and V can be easily determined by varying the direction of the relative velocity of the particle.Since V in ( 12) provides an upper bound on the velocity of any ejected particle, the physical meaning of this quantity is the gas outflow velocity.
In Figure 6, graphs of R 2 (t, 0, V ) corresponding to ejection at the perihelion of Phaethon's orbit, for two values of V = 1.23 km/s (red curve) and V = 1.76 km/s (blue curve), are superimposed on the graphs of S 2 (t).For V < 1.23 km/s, the values of R 2 are less than S 2 for any t, thus 1.23 km/s is a lower bound for V , assuming ejection at the perihelion.At speeds below 1.76 km/s, the stream's age is bounded above by the point t = 1600 years where the graphs of S 2 and its upper limit R 2 intersect.For higher speeds, R 2 increases and exceeds S 2 around t = 11 400 years -slightly to the left of the second minimum of deviation.In this case, the range of acceptable age values splits into two intervals: the first starting at t = 0 and the second in the vicinity of t = 11 400.
Thus, assuming the birth of the stream at the perihelion of Phaethon's orbit, the gas outflow velocity should have been at least 1.23 km/s.If we assume that this velocity did not exceed 1.76 km/s, the age of the stream is limited to a maximum of 1600 years.In general, age and ejection velocity estimates turn out to be related by the inequality ( 12) and computed values of S 2 (t): specifying one value gives bounds on the other.The region of permissible pairs of age and velocity values is highlighted in Figure 6.
If the ejection occurred outside the perihelion of the orbit, the estimates of velocities and age change.Calculations for the true anomaly value θ = 90 • , which corresponds to approximately twice the distance of the asteroid from the Sun compared to the perihelion, give a lower estimate of the gas outflow velocity equal to 1.06 km/s.If the velocity did not exceed 1.55 km/s, the age of the stream is estimated to be at most 1700 years.At higher speeds, as in the case of ejection at perihelion, the range of acceptable age values adds the vicinity of the second minimum of S 2 .

Robustness of results to sample variations
The age estimates of the Geminid stream and the gas ejection velocity discussed above are based on the 2020 meteoroid orbit sample.After simulating the orbit evolution, some of them were classified as strongly perturbed and discarded.A criterion for strong perturbation was the modulus of the five-year difference of one of the Keplerian elements falling into the upper 2% of the total distribution of such differences across the entire sample.The numerical values of critical differences are given in Section 3.2.This section addresses the question of how strongly the age and velocity values depend on the selectivity of the filter and the sample year.
Table 2 shows the values of τ 1 , T 2 , and V min calculated for different values of the filter selectivity α in the range from 10% to 0%.Recall that τ 1 is the age estimate obtained in Section 4.2, as the minimum of the dispersion S 5 , T 2 is the time interval between the minima of distances ϱ 5 and ϱ 2 between the average orbit of the sample and Phaethon's orbit, and V min is the lower estimate of the gas ejection velocity obtained assuming a cometary scenario of the stream's birth with Phaethon as the parent body.Additionally, the table includes the value V 1600 -the lower estimate of the gas outflow velocity for an age value of 1600 years.The column 'filtered out' indicates the percentage of orbits rejected by the filter.Age estimates in this and the following Table 3 are not rounded to centuries, as in the previous text, since in this case deviations, albeit small ones, from typical values are of interest.
Remarkably, the qualitative picture illustrated by the plots in Figures 4, 5, and 6 is preserved for all values of α, including zero, meaning no outlier filter.Quantitatively, we find the rows with α = 5%, 2%, 1% to be reliable.In the case of an element-wise threshold α = 10%, a quarter of the sample is filtered out, which is unacceptably high.In the unfiltered sample (α = 0%), the presence of significant outliers is visually evident in the plot of distances between individual orbits and the average.Non-robust statistics -mean and dispersion -are biased by orbits that have effectively left the stream due to perturbations.
α (%) τ 1 (years) T 2 (years) V min (km/s) V 1600 (km/s) filtered-out (%) Table 3 shows the spread of results by the sample year.All samples used in this table passed filtering with a selectivity of α = 2%.Compared to the previous table, columns S 2 (0) and S 5 (0) have been added, containing the dispersion values at present time.Samples from all five years give a minimum S 5 close to 1600 years.The minimum distances between the mean orbits of the samples and Phaethon's orbit fluctuate within 1200 to 1600 years, except for the ϱ 5 minimum for the orbits of the year 2019.This is the smallest sample, and the mentioned distance is close to constant over the interval of 1500 to 500 years.It is worth noting the low initial values of dispersions S 2 (0) and S 5 (0) for the sample of the year 2022 compared to other years.The corresponding lower estimate of the ejection velocity for this sample is also lower than the others because V min is determined by the smallest value of S 2 .

Conclusion
The minimum value of the dispersion S 5 for the modeled bundle of orbits is reached at the point τ 1 = 1600 years.With an accuracy of up to a century, this moment is the same for the samples of initial values of five consecutive years of observations.Starting from the τ 1 point, the deviation grows, changing not too rapidly compared to the previous time interval.
The birth of the Geminid stream at the moment τ 1 is a natural explanation for the minimization of S 5 .This hypothesis is supported by the result of Section 4.3, according to which the distance between the average orbit of the sample and Phaethon's orbit reaches a minimum over the interval T 2 = [1200,2000] years in the metrics ϱ 2 and ϱ 5 .
Assuming the formation of the Geminids as a result of the rapid destruction of a cometary nucleus, the remnant of which is Phaethon, the gas ejection velocity is estimated from below by the values υ 1 = 1.23 km/s for an ejection at perihelion and υ 2 = 1.06 km/s for an ejection at a point on the orbit with a true anomaly of 90 • .
The results of Section 4.4 establish the dependence between the lower limit of the gas ejection velocity and the permissible age values, presented in the diagram in Figure 6 for the case of gas ejection at perihelion.In particular, for velocities not exceeding 1.76 km/s and the ejection at perihelion, the corresponding age estimate is T 3 = [0, 1600] years.For the ejection at a point with θ = 90 • and a gas velocity of no more than 1.55 km/s, the age estimate is T 4 = [0, 1700] years.

Figure 1 :
Figure 1: Histograms of 99% absolute values of five-year differences of elements of the 2020 sample.The red dash marks the 98 percentile, that is the filter cutoff threshold

Figure 2 :
Figure 2: Dependence of elements of the mean orbit of the sample on time

Figure 3 :Figure 4 :
Figure 3: Dependence of standard deviations of angular orbital elements of the sample on time

Figure 5 :
Figure 5: Distances ϱ 2 and ϱ 5 between the mean orbit of the sample and Phaethon's orbit as a function of time.Black lines denote distance values, and colored lines represent moving averages with a two-sided window of ±250 years

Figure 6 :
Figure 6: The blue and red curves are the maximum possible values of S 2 , at a given gas outflow velocity equal to 1.76 and 1.23 km/s, respectively, at perihelion.The black curve is the actual S 2 values for the sample.The figure on the right shows the domain of permissible combinations of the ejection velocity and the age.

Table 1 :
Threshold differences of elements for the filters with selectivity 2% for samples of different years to merge datasets from different years to avoid biases in statistical characteristics due to varying observation conditions and instruments across different years.Let us note briefly, that qualitatively, the results are replicated in the other four datasets, and the corresponding quantitative differences are discussed in Section 4.5 of the article.

Table 2 :
Estimates of the age of Geminids and the gas outflow velocity for the 2020 sample with different filter selectivity α

Table 3 :
Estimates of the age of Geminids and the gas outflow velocity for samples from different years.Age estimates are in years, speeds are in km/s