Age of Geminids derived from the statistics of meteoroid orbits

ABSTRACT

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 non-linearly depends on the choice of orbital elements. In this work, using the construction of the Fréchet 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 analogues 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 yr in the past were conducted. By analysing 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. The estimate of the age of the stream lies in the interval from 1200 to 2400 yr, and the speed of gas outflow at perihelion should have been more than 1.2 km s⁻¹.

methods: data analysis, methods: statistical, celestial mechanics, meteorites, meteors, meteoroids, minor planets, asteroids: individual: (3200) Phaethon

1 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

$$\begin{eqnarray} \varrho _2(\boldsymbol{x}, \boldsymbol{y}) = \sqrt{(\boldsymbol{u}_x-\boldsymbol{u}_y)^2 + (\boldsymbol{v}_x-\boldsymbol{v}_y)^2}, \end{eqnarray}$$

(1)

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

$$\begin{eqnarray} \boldsymbol{u} = \frac{\boldsymbol{h}}{\varkappa },\quad \boldsymbol{v} = \boldsymbol{e}|\boldsymbol{u}|, \qquad \boldsymbol{h} = \boldsymbol{r} \times \boldsymbol{\dot{r}}, \quad \boldsymbol{e} = \frac{\boldsymbol{\dot{r}}\times \boldsymbol{h}}{\varkappa^2}-\frac{\boldsymbol{r}}{|\boldsymbol{r}|}, \end{eqnarray}$$

(2)

where ϰ² is the gravitational parameter, and r and |$\boldsymbol{\dot{r}}$| 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 non-trivially topologically structured.

In the context of meteoroid streams, the metric (1) is used both to search for parent bodies (Kokhirova et al. (2018); Sergienko, Sokolova & Kholshevnikov (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 ϱ₂.

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 outflow 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 yr 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.

2 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 non-linearity of the transformations between elemental systems.

This drawback of element-wise averaging of orbits was noted in the article Jopek, Rudawska & Pretka-Ziomek (2006). It was suggested to calculate the mean orbit as the best least-squares approximation of orbits in the family, parametrized 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 & 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

$$\begin{eqnarray} S^2(\boldsymbol{x}) = \frac{1}{n} \sum \limits_{k=1}^n D^2(\boldsymbol{x}_k, \boldsymbol{x}). \end{eqnarray}$$

(3)

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 & Molchanov (2005)), and in the case of a Riemannian manifold, uniqueness conditions are derived (Bhattacharya, Patrangenaru et al. (2003)), along with a generalization of the central limit theorem (Bhattacharya & Patrangenaru (2005)).

The element-wise mean orbit mentioned above falls within the definition of the Fréchet mean if we set

$$\begin{eqnarray} D^2(\boldsymbol{x}_1, \boldsymbol{x}_2) = \sum \limits_{i=1}^5\left(\backepsilon_i^1-\backepsilon_i^2\right)^2 , \end{eqnarray}$$

(4)

where |$\backepsilon$|₁, …, |$\backepsilon$|₅ are the elements of the orbit in the chosen system. The mean from Jopek et al. (2006) corresponds to the metric

$$\begin{eqnarray} D^2(\boldsymbol{x}_1, \boldsymbol{x}_2) = (\boldsymbol{h}_1-\boldsymbol{h}_2)^2 + (\boldsymbol{e}_1-\boldsymbol{e}_2)^2 + (E_1-E_2)^2. \end{eqnarray}$$

(5)

The value of the function (3) at the Fréchet mean, that is, the minimum over the entire orbit space |$\mathbb {H}$|⁠,

$$\begin{eqnarray} S^2 :=\min _{\boldsymbol {x}\in \boldsymbol{H}} S^{2}(\boldsymbol{x}) \end{eqnarray}$$

(6)

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 ϱ₂, it is possible to obtain exact formulas.

2.1 Mean with respect to the ϱ₂ metric

From definitions (3) and (1), it follows that to find the Fréchet mean in the ϱ₂ metric for a family of orbits defined by vectors (u_k, v_k), k = 1, …, n, one needs to minimize the quadratic form

$$\begin{eqnarray} S^2(\boldsymbol{u}, \boldsymbol{v}) = \frac{1}{n} \sum \limits _{k=1}^n\left((\boldsymbol{u}_k-\boldsymbol{u})^2 + (\boldsymbol{v}_k-\boldsymbol{v})^2\right), \end{eqnarray}$$

(7)

defined in |$\mathbb {R}^6$|⁠. Vectors u and v define a curvilinear orbit only if u ≠ 0 and

$$\begin{eqnarray} \boldsymbol{u} \boldsymbol{v} = 0. \end{eqnarray}$$

(8)

Thus, the form (7) should be minimized subject to the condition (8). The solution to this problem, as shown in Appendix  A, yields the following values for the vectors u and v of the mean orbit:

$$\begin{eqnarray} \boldsymbol{u} = \frac{1}{1-\mu ^2}\left(\bar{\boldsymbol{u}}-\mu \bar{\boldsymbol{v}}\right), \qquad \boldsymbol{v} = \frac{1}{1-\mu ^2}\left(\bar{\boldsymbol{v}}-\mu \bar{\boldsymbol{u}}\right), \end{eqnarray}$$

(9)

where

$$\begin{eqnarray} \bar{\boldsymbol{u}} = \frac{1}{n} \sum \limits _{k=1}^n \boldsymbol{u}_k, \quad \bar{\boldsymbol{v}} = \frac{1}{n} \sum \limits _{k=1}^n \boldsymbol{v}_k, \quad \mu = \frac{\left(|\bar{\boldsymbol{u}} + \bar{\boldsymbol{v}}| -|\bar{\boldsymbol{u}}-\bar{\boldsymbol{v}}|\right)^2}{4\bar{\boldsymbol{u}}\bar{\boldsymbol{v}}}. \end{eqnarray}$$

(10)

In exceptional cases where |$\bar{\boldsymbol{u}}$| is collinear with |$\bar{\boldsymbol{v}}$|⁠, μ² calculated by formulas (10) equals 1. In these cases (more details are discussed in Appendix  A), 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 (9),

$$\begin{eqnarray} \boldsymbol{u} \approx \bar{\boldsymbol{u}}-\mu \bar{\boldsymbol{v}}, \qquad \boldsymbol{v} \approx \bar{\boldsymbol{v}}-\mu \bar{\boldsymbol{u}} \end{eqnarray}$$

(11)

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₂ and express it in terms of the elements of the family:

$$\begin{eqnarray} S^2_{2}:=S^2(\boldsymbol{u}, \boldsymbol{v}) =\overline{\boldsymbol{u}^2} + \overline{\boldsymbol{v}^2}-\frac{1}{4}\left(|\bar{\boldsymbol{u}} + \bar{\boldsymbol{v}}| +|\bar{\boldsymbol{u}}-\bar{\boldsymbol{v}}|\right)^2, \end{eqnarray}$$

(12)

where

$$\begin{eqnarray} \quad \overline{\boldsymbol{u}^2} = \frac{1}{n}\sum \limits _{k=1}^n\boldsymbol{u}_k^2, \quad \overline{\boldsymbol{v}^2} = \frac{1}{n}\sum \limits _{k=1}^n\boldsymbol{v}_k^2. \end{eqnarray}$$

(13)

The relationship between the Keplerian orbit elements and the components of vectors u and v is determined by equalities (Kholshevnikov et al. (2016))

$$\begin{eqnarray} && \boldsymbol{u} = \sqrt{p} {\begin{pmatrix}\sin i\sin \Omega -\sin i\cos \Omega \\\cos i \end{pmatrix}},\\&& \boldsymbol{v} = e\sqrt{p} {\begin{pmatrix}\cos \omega \cos \Omega -\cos i\sin \omega \sin \Omega \\\cos \omega \sin \Omega +\cos i\sin \omega \cos \Omega \\\sin i\sin \omega \end{pmatrix}}, \end{eqnarray}$$

(14)

where e, p, i, ω, Ω are the eccentricity, semilatus rectum, inclination, argument of pericentre, 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, ω, Ω.

2.2 Mean with respect to to the ϱ₅ metric

When studying families of orbits that have undergone long-term evolution, the metric ϱ₂ may be too subtle a tool. Under the influence of small perturbations, the argument of pericentre 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

$$\begin{eqnarray} \varrho _{5} = \min \limits _{\Omega , \omega }\varrho _2. \end{eqnarray}$$

(15)

The distance ϱ₅ between orbits x and y is equal to the minimum distance ϱ₂ between orbits |$\mathcal {X}$| and |$\mathcal {Y}$|⁠, all Keplerian elements of which, except for Ω and ω, coincide with the corresponding elements of x and y.

Obviously, ϱ₅(x, y) = 0 if x and y differ only in the argument of pericentre and longitude of the ascending node. Thus, the orbit space |$\mathbb {H}$| is divided into equivalence classes: for any two elements of the same class, ϱ₅ equals zero, and formula (15) defines the distance between the classes. In Kholshevnikov et al. (2016), it is shown that ϱ₅ satisfies the triangle inequality. That is, the described set of equivalence classes, equipped with the distance ϱ₅, forms a metric space.

From the results of Milanov (2018), it follows that this space isometrically embeds into Euclidean space |$\mathbb {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.

$$\begin{eqnarray} && p = \frac{1}{n^2}\left(\sum _{j = 1}^n p_j + 2\sum _{1\leqslant j \lt k\leqslant n} \sqrt{p_jp_k}\cos(i_j-i_k)\right), \\ && e = \frac{1}{n\sqrt{p}}\sum _{j = 1}^n e_j\sqrt{p_j},\\ &&\cos i = \frac{1}{n\sqrt{p}}\sum _{j = 1}^n \sqrt{p_j}\cos i_j, \end{eqnarray}$$

(16)

where p_j, e_j, i_j, j = 1, …, n are the semilatus 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 (16) 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 (6), is expressed in Keplerian elements as follows:

$$\begin{eqnarray} S^2_5 = \frac{1}{n}\sum _{j = 1}^n p_j\big(1+e_j^2\big)-p(1+e^2). \end{eqnarray}$$

(17)

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₂ and S₅ 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 (6), Chebyshev’s inequality holds in the form

$$\begin{eqnarray} \mathbb {P}\left(D(x, \mu) \geqslant kS\right) \leqslant \frac{1}{k^2}, \end{eqnarray}$$

(18)

where |$\mathbb {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₂ and S₅ on time. We built these dependencies for the Geminid stream, modelling the evolution of the orbits of several thousand meteoroids over the past 20 000 yr. 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₅, 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.

3.1 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 (Battams et al. (2022)). 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, Ip & Spohn (2019), Cukier & Szalay (2023)).

The Geminids are part of a larger complex of the 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, Avdyushev & Williams (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)).

Battams et al. (2020, 2022) provided the summary and analysis of observations of the dust trail following the orbit of (3200) Phaethon as seen by the Wide-field Imager for the Parker Solar Probe and estimated the mass of dust. The value obtained is higher than can be provided by the modern activity of the asteroid but lower than the current mass estimates for the Geminids. Cukier & Szalay (2023) modelled 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. The current rate of outgassing is quite low, 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, so the scenario of cometary formation cannot be unambiguously excluded. In this case, the thickness of the dust crust above the water–ice layer should not have exceeded 1 m about 1000 yr 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 & Obrubov (1984) provides estimates in the range of 1600–19 000 yr, depending on the particle density. Using numerical integration, Jones (1985) obtained a value of 6000 yr. In Gustafson (1989), an interval of 1000–2000 yr is given, and it is noted that with more significant influence of non-gravitational forces, the estimate will decrease. The value of 6000 yr is mentioned in Rendtel (2005). By comparing the stream model with observational data, Ryabova (1999) finds the best agreement for an age of 2000 yr.

3.2 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. (2020, 2021)), available on the GMN website (https://globalmeteornetwork.org/data/traj_summary_data/). The selection of orbits for modelling the evolution was done yearly, covering the period from 2019 to 2023, resulting in five separate data sets. Particles identified with the Geminids meteor shower by GMN’s internal algorithms were considered. In total, we investigated five data sets, 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.

The published elements of meteoroid orbits take into account the correction associated with the gravitational influence of the Earth during close approach preceding the moment of observation. The technique of such ‘Earth subtraction’ is described, for example, in Jenniskens et al. (2011) and consists of inversely integrating the perturbed meteoroid orbit over a short period of time. In order not to take into account the mentioned correction again, the position of meteoroids was recalculated by shifting along the Keplerian orbit to the Julian date ten days before the start of observation in the corresponding year. Then, the meteoroid’s motion was computed backward in time for 20 000 yr, with a time-step of 0.005 yr.

The integration of motion equations was performed using the numerical integrator collo (https://github.com/shvak/collo), developed by V.Sh. Shaidulin. The simulation took into account the gravitational field of the Sun, eight planets, the Moon, and Pluto. The initial data and gravitational parameters for all listed bodies were taken from the EPM2021 catalogue of planetary ephemeris by IAA RAS (https://iaaras.ru/en/dept/ephemeris/epm/). 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 five-year 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 modelling period were excluded from further analysis. An orbit was considered significantly perturbed if the absolute value of the change in one of its elements during a five-year period exceeded the upper 2 per cent quantile of all such changes for the respective data set. Table 1 presents the threshold values for the absolute differences in elements, and Fig. 1 shows histograms of the 99 per cent absolute values of differences for the 2020 data set, indicating the cut-off threshold. This filtering procedure excludes 5 to 6 per cent of the orbits from each year’s data set.

Figure 1.

Histograms of 99 per cent absolute values of five-year differences of elements of the 2020 sample. The red dash marks the 98 percentile, which is the filter cut-off threshold.

Table 1.

Threshold differences of elements for the filters with selectivity 2 per cent for samples of different years.

Year	Δa (au)	Δe	Δi°	ΔΩ°	Δω°
2019	0.0130	0.00108	0.304	0.739	0.762
2020	0.0142	0.00126	0.315	0.825	0.842
2021	0.0134	0.00122	0.323	0.796	0.825
2022	0.0140	0.00127	0.316	0.751	0.791
2023	0.0146	0.00132	0.327	0.775	0.793

The results described further were obtained for the 2020 data set. We opted not to merge data sets 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 data sets, and the corresponding quantitative differences are discussed in Section 4.5 of the article.

4 RESULTS

4.1 Mean orbit of the sample

The evolution of Keplerian elements of the mean orbit, calculated using formulas (9), is shown in Fig. 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 semimajor axis undergoes changes of no more than 0.05 au until approximately 15 000 yr, after which this parameter sharply increases. The argument of pericentre, and longitude of the ascending node change more significantly than other elements, completing almost a full revolution over the simulated time span.

Figure 2.

Dependence of elements of the mean orbit of the sample on time.

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 Fig. 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. The arcsine of the eccentricity of an ellipse has a clear geometric meaning: it is the angle at which a segment connecting the centre with one of the foci 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 ϱ₂ metric, we will also consider the mean orbit in the ϱ₅ metric and the associated dispersion value S₅ defined by formula (17).

Figure 3.

Dependence of standard deviations of angular orbital elements of the sample on time.

Recall that ϱ₅ 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 pericentre. The mean orbits with respect to the ϱ₂ and ϱ₅ metrics are close to each other for most of the simulation period. The distance ϱ₂ between them does not exceed 0.005 for the first 15 000 yr and then sharply increases.

4.2 Dispersion of the orbits sample

Let’s examine the evolution of quantities S₂ and S₅ over time. Defined by the general formula (6), these values represent analogues of sample variance on the manifold of orbits equipped with metrics ϱ₂ and ϱ₅. Fig. 4 depict the evolution of S₂ and S₅ over the entire simulation period and during the first 4000 yr, which is of particular interest. Additionally, for the 2020 sample, orbit modelling was extended 10 000 yr into the future. The plots of S₂ and S₅ over the interval −20 000 to 10 000 are shown in the third panel of the figure. The minimum of S₅ occurs at the 1600 yr mark; hereafter, we will round time values to the nearest century. At this time, the orbits of the simulated stream were most densely concentrated around the mean in the ϱ₅ metric. According to our hypothesis, the time τ₁ = 1600 yr is the best candidate for the age of the Geminids, representing the global minimum of S₅ closest to the present.

Figure 4.

Dependence of dispersions S₂ and S₅ of sample orbits on time for different time scale fragments.

The intervals of monotonicity for the quantity S₂ closely resemble the corresponding pattern for S₅. However, the minimum of S₂ is located later – around the 300-yr 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 Fig. 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.

4.3 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 modelled over the same time span. Plots of distances ϱ₂ and ϱ₅ between the mean stream orbit and the asteroid’s orbit are presented in Fig. 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 coloured curves depict the moving averages with a window of ±250 yr.

Figure 5.

Distances ϱ₂ and ϱ₅ between the mean orbit of the sample and Phaethon’s orbit as a function of time. Black lines denote distance values, and coloured lines represent moving averages with a two-sided window of ±250 yr.

The global minima of the averaged dependencies are located at 1700 yr for ϱ₂ and 1200 yr for ϱ₅. Thus, according to the modelling results, during the interval T₂ = [1200, 1700] yr ago, the considered stream of orbits had the smallest size in terms of the root mean square deviation S₅, and its mean orbit was closest to the orbit of the presumed parent body.

The distance ϱ₅ also exhibits a deep local minimum at 10 000 yr, but the nearest minima of S₂ and S₅ differ by more than 2000 yr.

4.4 Dispersion of the orbits sample and the gas outflow velocity

It is interesting to compare the obtained values of S₂ 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 gas outflow speed 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 ϱ₂ 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 0.500 to 0.83 km s⁻¹ for dust particles and 0.945 to 1.368 km 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 |$\dot{\boldsymbol {r}}(t, \theta)$|⁠. We introduce the quantity R₂ = R₂(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 |$\dot{\boldsymbol {r}}(t, \theta)$| by no more than V in magnitude. From the definition (6) of dispersion, it follows that at the moment of ejection, for any sample of meteoroid orbits X of size n:

$$\begin{eqnarray} S_2 &=& \left(\frac{1}{n}\min _{m \in \mathbb {H}}\sum _{x \in X} \varrho _2^2(x, m)\right)^{1/2} \leqslant \left(\frac{1}{n}\sum _{x \in X} \varrho _2^2(x, x_p)\right)^{1/2}\\ &\leqslant& R_2(t, \theta , V), \end{eqnarray}$$

(19)

assuming that the particle velocities do not exceed V.

The value of R₂ for given t, θ, and V can be easily determined by varying the direction of the relative velocity of the particle. Since V in (19) provides an upper bound on the velocity of any ejected particle, the physical meaning of this quantity is the gas outflow velocity.

In Fig. 6, graphs of R₂(t, 0, V) corresponding to ejection at the perihelion of Phaethon’s orbit, for two values of V = 1.23 (red curve) and 1.76 km s⁻¹ (blue curve), are superimposed on the graphs of S₂(t). For V < 1.23 km s⁻¹, the values of R₂ are less than S₂ 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 yr, where the graphs of S₂ and its upper limit R₂ intersect. For higher speeds, R₂ increases and exceeds S₂ around |$t = 11\, 400$| yr – 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$|⁠.

Figure 6.

The blue and red curves are the maximum possible values of S₂, 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₂ values for the sample. The figure on the right shows the domain of permissible combinations of the gas outflow velocity and the age.

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 yr. In general, age and gas outflow velocity estimates turn out to be related by the inequality (19) and computed values of S₂(t): specifying one value gives bounds on the other. The region of permissible pairs of age and velocity values is highlighted in Fig. 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 yr. At higher speeds, as in the case of ejection at perihelion, the range of acceptable age values add the vicinity of the second minimum of S₂.

The obtained estimates of the gas outflow velocity are high but do not exceed reasonable limits. In the works Kasuga & Masiero (2022) and Zhang et al. (2023), the sodium outflow velocity from Phaethon is estimated at 0.7–0.8 and 1 km s⁻¹, respectively. Ryabova (2013) gives an interval of about 0.9–1.4 km s⁻¹. Let us recall that our estimates are obtained under the assumption of rapid formation of the Geminids due to the intense destruction of the cometary nucleus.

4.5 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 per cent 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 τ₁, T₂, and V_min calculated for different values of the filter selectivity α in the range from 10 to 0 per cent. Recall that τ₁ is the age estimate obtained in Section 4.2, as the minimum of the dispersion S₅, T₂ is the time interval between the minima of distances ϱ₅ and ϱ₂ 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₁₆₀₀ – the lower estimate of the gas outflow velocity for an age value of 1600 yr. 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.

Table 2.

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

α (per cent)	τ1 (years)	T2 (years)	Vmin (km s−1)	V1600 (km s−1)	filtered-out (per cent)
10	1730	1680–2330	1.22	1.56	25
5	1660	1370–1950	1.23	1.66	13
2	1600	1210–1680	1.23	1.78	6
1	1560	1205–1530	1.23	1.84	3
0	1370	780–1375	1.23	1.98	0

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⁻¹.

Year	τ1	T2	Vmin	V1600	filtered-out (per cent)	N	S2(0)	S5(0)
2019	1575	775–1845	1.24	1.71	6	2116	0.0291	0.0215
2020	1600	1210–1680	1.23	1.78	6	5959	0.0281	0.0225
2021	1590	1465–1985	1.23	1.78	5	9968	0.0281	0.0215
2022	1625	1520–1815	1.0	1.51	6	15 795	0.0235	0.0186
2023	1615	1605–1955	1.18	1.67	6	19 653	0.0266	0.0202

Remarkably, the qualitative picture illustrated by the plots in Figs 4, 5, and 6 is preserved for all values of α, including zero, meaning no outlier filter. Quantitatively, we find the rows with |$\alpha = 5, 2, 1\ \hbox{per cent}$| to be reliable. In the case of an element-wise threshold, |$\alpha =10\ \hbox{per cent}$|⁠, a quarter of the sample is filtered out, which is unacceptably high. In the unfiltered sample (⁠|$\alpha =0\ \hbox{per cent}$|⁠), 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.

Table 3 shows the spread of results by the sample year. All samples used in this table passed filtering with a selectivity of |$\alpha =2\ \hbox{per cent}$|⁠. Compared to the previous table, columns S₂(0) and S₅(0) have been added, containing the dispersion values at present time. Samples from all five years give a minimum S₅ close to 1600 yr. The minimum distances between the mean orbits of the samples and Phaethon’s orbit fluctuate within 1200 to 1600 yr, except for the ϱ₅ 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 yr. It is worth noting the low initial values of dispersions S₂(0) and S₅(0) for the sample of the year 2022 compared to other years. The corresponding lower estimate of the gas outflow velocity for this sample is also lower than the others because V_min is determined by the smallest value of S₂.

5 CONCLUSION

The minimum value of the dispersion S₅ for the modelled bundle of orbits is reached at the point τ₁ = 1600 yr. 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 τ₁ point, the deviation grows, changing not too rapidly compared to the previous time interval.

The birth of the Geminid stream at the moment τ₁ is a natural explanation for the minimization of S₅. 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₂ = [1200, 2000] yr in the metrics ϱ₂ and ϱ₅.

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 outflow velocity is estimated from below by the values υ₁ = 1.23 km s⁻¹ for an ejection at perihelion and υ₂ = 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 Fig. 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₃ = [0, 1600] yr. 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₄ = [0, 1700] yr.

The samples of meteoroid orbits whose evolution simulations produce the above age and gas outflow speed estimates represent only a small portion of the stream’s orbits, since all of these trajectories encounter the Earth. How much influence can the selection effect have on the results obtained? The method for estimating gas velocity described in Section 4.4 is robust to this effect. Indeed, the upper bound for the variance R₂ does not depend on the sample, and the inequality (19) is valid for any set of orbits. Stream age estimates obtained in Sections 4.2 and 4.3 are more vulnerable to selection. However, these two estimates, obtained by different methods on the same data, are consistent with each other. The characteristic size of the sample, described by the variance S₅, was smallest at the same time when the average orbit of this sample was closest to the orbit of Phaethon. The consistency of the estimates speaks in favour of their unbiasedness, which is further confirmed by the results of 4.4, which require a very high gas velocity for ages greater than 1600 yr.

ACKNOWLEDGEMENTS

We are grateful to the anonymous referee for the useful comments, which improved the article.

DATA AVAILABILITY

The data underlying this article are available at https://disk.yandex.ru/d/HLWHRZzcQwLVwQ. These data sets were derived from sources in the public domain: Global Meteor Network website (https://globalmeteornetwork.org/data/traj_summary_data/) and the online catalogue of planetary ephemeris EPM2021 IAA RAS (https://iaaras.ru/en/dept/ephemeris/epm/). The numerical integrator program code is availiable at https://github.com/shvak/collo.

APPENDIX: Calculation of vector elements of the mean orbit in the metric ϱ₂

Minimizing the expression (7) under the constraint (8) reduces to minimizing the Lagrangian function L(u, v, λ) = S(u, v) + λuv. The critical points of L satisfy the system of equations

$$\begin{eqnarray} \boldsymbol{u} + \mu \boldsymbol{v} = \bar{\boldsymbol{u}}, \qquad \boldsymbol{v} + \mu \boldsymbol{u} = \bar{\boldsymbol{v}}, \qquad \boldsymbol{u} \boldsymbol{v} = 0, \end{eqnarray}$$

(A1)

where

$$\begin{eqnarray} \bar{\boldsymbol{u}} = \frac{1}{n} \sum \limits _{k=1}^n \boldsymbol{u}_k, \qquad \bar{\boldsymbol{v}} = \frac{1}{n} \sum \limits _{k=1}^n \boldsymbol{v}_k, \qquad \mu = \frac{\lambda }{2n}. \end{eqnarray}$$

(A2)

The solutions u, v are given by the formulas

$$\begin{eqnarray} \boldsymbol{u} = \frac{1}{1-\mu ^2}\left(\bar{\boldsymbol{u}}-\mu \bar{\boldsymbol{v}}\right), \qquad \boldsymbol{v} = \frac{1}{1-\mu ^2}\left(\bar{\boldsymbol{v}}-\mu \bar{\boldsymbol{u}}\right), \end{eqnarray}$$

(A3)

and the parameter μ is zero if |$\bar{\boldsymbol{u}}\bar{\boldsymbol{v}} = 0$|⁠, or it is found from the equation

$$\begin{eqnarray} \mu ^2-\mu \nu + 1 = 0, \qquad \nu = \frac{\bar{\boldsymbol{u}}^2 + \bar{\boldsymbol{v}}^2}{\bar{\boldsymbol{u}}\bar{\boldsymbol{v}}}. \end{eqnarray}$$

(A4)

The second differential of L

$$\begin{eqnarray} {\rm d}^2L = 2(({\rm d}\boldsymbol{u})^2 + ({\rm d}\boldsymbol{v})^2 + 2\mu\,{\rm d}\boldsymbol{u}\,{\rm d}\boldsymbol{v}) \end{eqnarray}$$

(A5)

is transformed into the canonical form by changing variables du = dx + dy, dv = dx − dy

$$\begin{eqnarray} {\rm d}^2L = 4((1+\mu)({\rm d}\boldsymbol {x})^2 + (1-\mu)({\rm d}\boldsymbol {y})^2). \end{eqnarray}$$

(A6)

Hence, the minima of S are calculated using the formulas (A3), (A4), under the condition |μ| ≤ 1. Such a restriction is satisfied by exactly one root of the equation (A4):

$$\begin{eqnarray} \mu = \frac{1}{2}\left(\nu -\operatorname{sgn}{\nu } \sqrt{\nu ^2-4}\right) = \frac{\left(|\bar{\boldsymbol{u}} + \bar{\boldsymbol{v}}| -|\bar{\boldsymbol{u}}-\bar{\boldsymbol{v}}|\right)^2}{4\bar{\boldsymbol{u}}\bar{\boldsymbol{v}}}. \end{eqnarray}$$

(A7)

The formulas (A3) do not work in the case |$\bar{\boldsymbol{u}} = k\bar{\boldsymbol{v}}$| for a real k, |k| ≤ 1. If |k| < 1, then expressions (A7) and (A3) give μ = k and u = 0. Recall that straight-line orbits are not part of the space |$\mathbb {H}$| and, therefore, in this case, the problem has no solution.

If k = ±1, then μ = ±1 and the values of u, v are undefined. Denoting |$\boldsymbol{w} = \bar{\boldsymbol{u}} = k\bar{\boldsymbol{v}}$|⁠, and using the orthogonality condition (8), express S(u, v) as

$$\begin{eqnarray} && S(\boldsymbol{u}, \boldsymbol{v}) = \overline{\boldsymbol{u}^2} + \overline{\boldsymbol{v}^2}-\boldsymbol{w}^2 + (\boldsymbol{u} + k\boldsymbol{v}-\boldsymbol{w})^2,\\ &&\overline{\boldsymbol{u}^2} = \frac{1}{n}\sum \limits _{k=1}^n\boldsymbol{u}_k^2, \quad \overline{\boldsymbol{v}^2} = \frac{1}{n}\sum \limits_{k=1}^n\boldsymbol{v}_k^2. \end{eqnarray}$$

(A8)

Any pair of orthogonal vectors u, v that satisfies the equation u + kv = w minimizes S in this case.

Thus, except for the special cases described above, the Fréchet mean in the space of curved Keplerian orbits |$\mathbb {H}$| with the metric ϱ₂ exists, is unique, and is computed using the formulas (A3).