© The Authors 2024

1 Introduction

Stream identification and differentiation from the sporadic background are crucial components of meteor science. The International Astronomical Union (IAU) actively updates a detailed registry¹ of meteor showers and their originating celestial bodies. This registry is frequently augmented by findings from new meteor or fireball observation surveys, highlighting the evolving landscape of this scientific domain (Jenniskens et al. 2009; Kornoš et al. 2014; Rudawska et al. 2015; Jenniskens et al. 2016). Identifying these streams relies on employing a dissimilarity measure known as the D-function. Initially introduced by Southworth & Hawkins (1963), the D_SH is derived by comparing two groups of orbital elements. This measure, which increases as the similarity decreases, is often termed a “dissimilarity parameter.” Several adaptations have been introduced, each slightly changing the degree of influence of different elements to the D value (Drummond 1981; Steel et al. 1991; Jopek 1993; Jopek et al. 2008; Jenniskens et al. 2009).

Despite its extensive application, the efficacy of the D values or orbital similarity in general has been contested, particularly its dependency on sample size and the risk of falsely identifying random groupings as actual meteor showers (Pauls & Gladman 2005; Koten et al. 2014; Egal et al. 2017; Vida et al. 2018; Shober & Vaubaillon 2024). With the surge in accessible meteor orbit data from sources like video meteor networks, the demand for more reliable techniques to verify links between meteor showers and their parent bodies has intensified. This requirement is amplified by recent acknowledgments of the higher-than-anticipated uncertainty surrounding meteoroid orbital elements (Egal et al. 2017; Vida et al. 2018; Shober et al. 2023). To address these issues, new statistically robust methods have been proposed, such as the wavelet transform (Galligan & Baggaley 2002), application of the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) algorithm to identify clusters in meteor trajectory data (Moorhead 2016; Sugar et al. 2017) or the estimation of meteor shower false positive detections leveraging Kernel Density Estimation (KDE) (Shober & Vaubaillon 2024).

In spite of the ongoing challenges, at least 110 meteoroid streams² have been confirmed and solidly associated with primarily dust ejection from cometary activity (Jopek & Kaňuchová 2017; Jenniskens et al. 2020). Streams of larger meteoroids and asteroids have been the subject of many studies (Fu et al. 2005; Schunová et al. 2012; Jopek 2020); however, no significant evidence supports the associations between the meteoroids and NEOs. Recently, likely clusters within the NEO population have been identified by Jopek (2020); however, no considerable statistical evidence has supported associations between fireball data and such NEOs. If a macroscopic meteoroid stream existed, it could pose a consistent impact threat, and thus, the possibility has been continuously explored. Pauls & Gladman (2005) was the first comprehensive study aimed to investigate whether large meteoroid streams, proposed to be responsible for the similarity of a handful of meteorite falls, could remain coherent over long enough periods to be found in the fireball and meteorite fall datasets. Initially, they explored the statistical likelihood of orbital similarity within known meteorite-dropping fireballs, comparing these occurrences against randomly generated samples of orbits. Pauls & Gladman (2005) performed detailed numerical simulations on three well-documented cases of meteorite falls (Innisfree, Peekskill, and Příbram) that had been suggested to belong to meteoroid streams due to their orbital similarities with other falls or fireballs.

The objective of Pauls & Gladman (2005) was to address the three arguments used for the existence of macroscopic meteorite-dropping meteoroid streams: (1) orbital similarity (Příbram-Neuschwanstein and Innisfree-Ridgedale) (Halliday 1987; Spurnỳ et al. 2003; Halliday et al. 1990), (2) the clustering of fireball data (McCrosky et al. 1979; Halliday et al. 1990), and (3) meteorite falls which had the same day-of-fall and similar chemical trace elements (Lipschutz et al. 1997). To address these arguments, Pauls & Gladman (2005) determined how long these assumed streams could maintain their coherence before becoming indistinguishable from the sporadic background. They showed that if these three meteorites were associated with mete-oroid streams formed in near-Earth space, they would become decoherent on relatively short timescales of 10⁴ to 10⁵ years. This rapid dispersion would require an extremely recent breakup of the parent body, much more recent than the cosmic ray exposure ages of the meteorites involved and much more frequently than near-Earth asteroid impact rates (Bottke et al. 1994). Furthermore, the statistical evidence for clustering similar orbits was demonstrated to occur just as frequently in a random sample of orbits (Pauls & Gladman 2005).

Since these results nearly two decades ago, many further studies have tried to use the orbital similarity of meteorite falls to find parent bodies or streams. This warrants a renewed analysis of the statistical significance of clustering in the meteorite fall and fireball datasets. Since 2005, there has been an exponential growth in the number of meteor and fireball observation networks, significantly augmenting our knowledge and ability to detect smaller streams. This is exemplified by the fact that in 2005, there were 10 known meteorite falls with precise orbits, whereas in 2024, there are at least 55 in total, with 18 alone recovered in the last five years. This is primarily due to several new networks operating around the globe dedicated to recovering meteorites (Bland et al. 2012; Devillepoix et al. 2020; Colas et al. 2020; Borovička et al. 2022b). Given this increased capacity and the recent use of orbital similarity to make grand conclusions about source regions or parent bodies, we must re-address this question of macroscopic meteoroid streams.

2 Methods

To address the utility of orbital similarity for meteorite-dropping fireballs, near-Earth objects, and the existence of large meteoroid streams, we have used three primary methods: (1) a false-positive estimation and orbital clustering significance analysis using D_{S H} (Southworth & Hawkins 1963), D′ (Drummond 1980), and D_H (Jopek 1993), (2) Lyapunov characteristic lifetime mapping of a-e-l space, and (3) decoherence lifetime estimation of 300 fictitious streams.

2.1 Data collection and reduction

We examine several datasets in this study, including (1) 50 meteorite falls with orbits, (2) 616 potential 1 g meteorite-dropping fireballs collected by the FRIPON, EFN, and GFO networks (350 of which are probably 50 g droppers), (3) 3290 FRIPON fireballs, (4) 824 EFN fireballs taken from (Borovička et al. 2022b), (5) 35 012 telescopically observed NEOs from the JPL HORIZONS³ database, and (6) 310 impacts detected by USG sensors. The reliability of the USG sensor orbits has been demonstrated to be inconsistent and certainly not ideal for finding streams (Devillepoix et al. 2019); however, the data was included as recent studies have relied on the dataset using D-discriminants.

2.1.1 Fireball observations

The fireball data used in this work are derived from several key fireball observation networks, each contributing unique observations crucial for this analysis. The first of these is the Fireball Recovery and InterPlanetary Observation Network⁴ (FRIPON), which is an expansive network primarily covering Europe, equipped with over 220 all-sky cameras covering 15 countries across the globe. It focuses on both detecting fireballs and recovering meteorites with its CCD video observatories. FRIPON’s system is particularly designed to capture fireballs brighter than magnitude zero, providing detailed data for trajectory and orbit determination for meteorite-dropping fireballs along with some smaller, largely shower components (Colas et al. 2020).

We also utilized data collected by the Global Fireball Observatory⁵ (GFO). The GFO is an expansive collaboration that leverages the extensive infrastructure developed by the Desert Fireball Network (DFN) (Devillepoix et al. 2020). The partnership comprises ten partner networks, supported by 18 collaborating institutions spread over nine countries globally. There are over 100 observatories spread across the networks, with 50 covering over a third of the Australian landmass. The weather-hardened digital fireball observatories feature high-resolution DSLR cameras, all-sky fisheye lenses, and GNSS-synchronized liquid crystal shutters for precise timing. Facilitating the recovery of 15 meteorites since the DFN’s digital upgrade in 2013–2015, the collaboration has recovered approximately 30% of all meteorites recovered with known orbits to date (King et al. 2022; Shober et al. 2022; Devillepoix et al. 2022; Anderson et al. 2022). For trajectory analysis, the GFO employs a straight-line least squares (SLLS) method combined with an extended Kalman smoother to ensure accurate velocity estimates (Borovička 1990; Sansom et al. 2015). This sophisticated approach integrates observational uncertainties to manage and propagate errors effectively. Such methodologies allow for the precise determination of pre-entry orbits by numerically integrating meteoroids’ states beyond Earth’s gravitational influence, taking into account all significant perturbative effects (Jansen-Sturgeon et al. 2019; Shober et al. 2019, 2020).

Lastly, the European Fireball Network (EFN) is the longest-running photographic fireball network, with the first observations taking place in 1963. It has a total coverage area of about one million square kilometers. A history spanning several decades, the EFN team has successfully facilitated the recovery of at least 13 meteorites (Spurnỳ et al. 2003, 2020; Spurny et al. 2013; Spurný et al. 2017). The network primarily operates digital autonomous fireball bbservers (DAFOs), which are weather-proof and fully autonomous systems designed to continuously monitor the entire sky during clear weather conditions. Each DAFO utilizes DSLR cameras fitted with fisheye lenses, similar to those used by the DFN, capturing high-resolution images of meteor events. These cameras are configured to detect meteors with an absolute magnitude brighter than −2 and can perform detailed radiometric measurements for fireballs exceeding −4 in magnitude (Borovička et al. 2022b). EFN’s detection system is adept at identifying meteoroids of mass greater than 5 grams and capable of capturing high-velocity meteoroids as small as 0.1 grams. Spread across central Europe, the network’s 26 stations cover an area of approximately one million km², providing extensive observational data.

The dataset from EFN utilized in this study has been detailed in Borovička et al. (2022b) and Borovička et al. (2022a). The atmospheric trajectory of each fireball is deduced using the straight line least squares (SLLS) method, which assumes a linear trajectory through space (Borovicka 1990). Unlike the DFN, EFN employs a different approach for velocity calculation, using time data projected onto the computed SLLS trajectory and fitting it with a physical model that accounts for atmospheric drag and ablation (Pecina & Ceplecha 1983). This model includes parameters for preatmospheric velocity, ablation coefficient, and mass-related properties, refined using atmospheric density models like CIRA72 or NRLMSISE-00. Manual adjustments or alternative models are used when significant deceleration occurs. Outliers were evaluated or removed to ensure data accuracy, and any systematic discrepancies between cameras were resolved before finalizing the heliocentric orbits, which are calculated with a modified version of the method described in Ceplecha (1987), accounting for Earth’s rotation and gravitational influences.

The 350 “meteorite-dropping” fireballs observed by the FRIPON, EFN, and GFO projects were identified primarily using the α − β methodology of Sansom et al. (2019). The α − β coefficients correspond to the ballistic coefficient and mass-loss parameter, respectively. These parameters, introduced and developed in Gritsevich (2007) and Lyytinen & Gritsevich (2016), can be estimated for any impact event with some degree of deceleration. These parameters provide a simple yet effective way to characterize the probability of a likely fall without needing a multivariate solution or supercomputing resources (Sansom et al. 2019). The 616 potential dropping events were identified within the FRIPON and GFO datasets using a variation of Equations (7) and (8) in Sansom et al. (2019), using an end mass of 1 g (616 events) and 50 g (350 events). The EFN meteorite dropping subset was taken from Borovička et al. (2022b), using only events with an end mass above 1 g or 50 g according to their ablation and fragmentation model.

2.1.2 NEOdata

The NEO data utilized in this study was taken from the Jet Propulsion Laboratory (JPL) Horizons database, which provides comprehensive ephemeris computations essential for the precise location of celestial objects over time. Accessible via NASA’s Solar System Dynamics Group, the Horizons system offers a robust, highly reliable tool for generating accurate solar, lunar, and planetary ephemerides. Utilized extensively in the astronomical community, this service facilitates detailed orbital analysis by delivering key parameters such as positions, velocities, and magnitudes across various time points. These data are crucial for high-precision object-tracking, making Horizons an invaluable resource for observational planning and ongoing research in celestial mechanics.

2.1.3 USG sensors detections

The Center for Near-Earth Object Studies (CNEOS), operated by NASA’s JPL, also maintains a fireball database⁶ that catalogs detailed observations of fireballs produced by meteoric phenomena in Earth’s atmosphere. This database is compiled from data collected by U.S. Government sensors, offering vital information on the trajectory, speed, penetration depth in the atmosphere, and estimated energy release of these events. The CNEOS fireball database is a critical tool for understanding the population of large meteoroids in terms of impact energy distribution, providing global coverage for detecting large impactors. This database is very useful for increasing the statistical significance of our understanding of the subset of the largest annual impactors. In total, 310 CNEOS impact events to date have velocity and peak brightness latitude, longitude, and altitude information. Unfortunately, this valuable data is intentionally degraded to maintain the secrecy around the sensitive nature of the US military’s monitoring capabilities. Thus, the data released is very unreliable for orbit determinations, with velocities varying up to over 20% and the radiant, in some cases, being 90° off when compared to more precise ground-based fireball observations (Devillepoix et al. 2019). Therefore, using this data with orbital similarity functions is quite unideal and ill-advised. Despite this, we included this dataset given its extensive use in several studies. With that being said, we deliberately chose to not try to estimate the large uncertainties and to test the statistical significance of the CNEOS fireballs even if the errors were minuscule.

2.2 Statistical significance of orbital similarity

To assess the statistical significance of orbital similarities, we employed various D-functions to measure the dissimilarity between the orbital elements of celestial bodies. The functions applied included the D_{S H}, D′, and D_H, which consider different combinations of orbital parameters such as perihelion distance (q), eccentricity (e), inclination (ι), argument of perihelion (ω), and longitude of ascending node (Ω).

The analysis involved two aspects: (1) using the methodology of Shober & Vaubaillon (2024) to estimate the false positive clustering rate by fitting a KDE to the concerned populations, and (2) calculating the D values for all pair combinations within the concerned dataset and examining the distributions.

Using the KernelDensity class from scikit-learn, as done in Shober & Vaubaillon (2024), we estimated the false positive association rate and the statistical significance of alleged meteorite-dropping streams. KDE is a non-parametric technique that estimates the probability density function of a random variable from its data points. This is achieved by placing a kernel, often a Gaussian function, on each data point and aggregating these kernels to form a smooth representation of the underlying probability density function. This technique is adept at modeling unknown distributions and handling the multimodality frequently observed in sparse datasets (Silverman 2018). Moreover, KDE is recognized for its adaptability in accurately estimating densities across a range of shapes, assuming the smoothing parameter is correctly chosen (Seaman & Powell 1996).

We used KDE to estimate the orbital distribution of the NEO population, the USG sensor impact population, and the 50 meteorite falls with orbits. This involved generating synthetic datasets from these KDEs, equal in size to the original datasets, and evaluating the number of pairs below specific D-value thresholds to estimate false positives. This approach helps us understand if observed pairings in the data could be explained by random chance rather than physical association. A Monte Carlo simulation repeated this process 500 times for the meteorite-NEO and USG sensors-NEO pairs; from this, the standard deviations of the sample means for the numbers of random associations could be calculated to understand the significance of the number of actual associations between the populations. Please refer to Shober & Vaubaillon (2024) for a more detailed explanation of the KDE false-positive methodology.

Additionally, we calculated all possible D values for every pair combination within a single dataset. This extensive analysis allowed us to evaluate whether the D-value distribution followed expected patterns from random chance or exhibited clusters indicative of genuine associations. If a significant, actively occurring process was generating streams in near-Earth space, there should be an excess of small D values within the distributions. The total number of unique combinations within an orbital dataset follows $C=\frac{n(n-1)}{2}$, where n is the size of the dataset. In total, within this study, we have calculated the orbital similarity D values for over 10¹¹ unique orbital pairs between the meteorite-fireball-NEO comparisons and the Monte Carlo procedure using KDE to estimate the statistical significance of the similarities.

2.3 Lyapunov characteristic lifetimes

The Lyapunov characteristic exponent (LCE) quantitatively measures chaos in a dynamical system. It characterizes the rate at which trajectories initially close in phase space diverge over time, indicating the system’s sensitive dependence on initial conditions. This chaos is characterized by an exponential divergence. To calculate the LCE, we employ the open-source N-body simulation software, Rebound, specifically tailored for astronomical systems (Rein & Liu 2012). The solar system model included the Sun and the eight major planets with initial conditions sourced from JPL Horizons⁷. The ‘whfast’ integrator was utilized for its efficiency with near-Keplerian orbits, with a time step set to 0.01 years.

Orbits were non-uniformly generated across a parameter space defined by semi-major axes derived from the Tisserand parameter relative to Jupiter, with semi-major axis values ranging from 0.5 to 3.5 au, eccentricities ranging from 0 to 0.95, and an inclination of 0, 10, or 20 degrees. The non-uniform distribution of particles was chosen to concentrate the particle density in chaotic regions and optimize the resolving process. In total, 289640 particles at 0° inclination, 287 352 at 10° inclination, and 353 470 particles at 20° inclination were integrated for the map constructions (Sect. 3.3). A lower inclination range was chosen as the asteroidal population is nearly all within this range of values. The integration time for each simulation was set to 20 000 years, consistent with other studies (Tancredi 1995, 1998; Shober et al. 2024b), and the Lyapunov time was calculated as the inverse of the Lyapunov exponent.

The Rebound library integrates the dynamical equations of motion alongside the variational equations, following the methodology described by Cincotta & Simó (2000). These variational equations describe how infinitesimal perturbations to the system’s state vector (position and velocity) evolve over time. A description of this can be found within the documentation⁸ of Rebound. The integration process involves the following steps:

1. Numerically solving the N-body system’s equations of motion using an integrator like WHFast, IAS15, or SEI.

2. Simultaneously integrating the variational equations to track how perturbations in initial conditions evolve.

3. Applying orthonormalization regularly (e.g., Gram-Schmidt process) to the deviation vectors to maintain numerical stability and prevent scaling issues.

4. Estimating the Lyapunov exponent from the average exponential rate of divergence of these orthonormalized vectors, normalized by the integration time.

This method enables efficient and accurate estimation of the Lyapunov characteristic lifetime, which is vital for understanding the stability and evolution of meteoroid streams under gravitational influences. For a detailed description of the integration methods and theoretical background, refer to the work by Cincotta & Simó (2000) and Rein & Liu (2012).

2.4 Meteoroid stream decoherence lifetimes

Meteoroid stream decoherence lifetimes were determined by simulating the dynamical evolution of stream-like groups of particles under the solar system’s gravitational influences. The simulations were also performed using the N-body integration software REBOUND, configured with the solar system’s major planets as perturbing bodies and using the IAS15 integrator with an adaptive timestep (Rein & Liu 2012; Rein & Spiegel 2015). Streams were modeled as packets of meteoroids, each originating from similar orbital elements but including small initial variations to simulate natural dispersion. The initial positions of the 100 particles in each packet were identical, but an initial uniform velocity variation between 0–5 m s⁻¹ was added in each x–y–z direction. The ejection point true anomaly on the orbit was also chosen from a uniformly random distribution. For each packet, the primary meteoroid defined the base orbit, with subsequent meteoroids perturbed slightly in velocity space to represent typical ejection velocities from a parent body. The semi-major axis and eccentricity space were mapped out with a grid layout, allowing the generation of 100 fictitious streams per inclination angle (0°, 10°, and 20°). The simulation data was used to construct decoherence maps, illustrating the decoherence lifetimes as functions of semi-major axis and eccentricity for different inclinations.

The dynamical evolution of these particles was tracked over a simulation time of 5 million years, as the decoherence lifetimes of streams tend to be on the scale of 10⁴–10⁵ years (Pauls & Gladman 2005). The spread of orbital elements over time was monitored, focusing on the size of the largest coherent cluster within each packet. The DBSCAN (density-based spatial clustering of applications with noise) algorithm was employed to identify and track the sizes of these clusters as it has been previously demonstrated as an effective way to identify clusters in meteor datasets (Moorhead 2016; Sugar et al. 2017). The parameters for the DBSCAN algorithm included:

• Epsilon (ϵ): The maximum distance between two points to be considered as neighbors. This value varied based on the specific D-function used and was based on the results of the NEO similarity significance analysis:

  • D_SH: ϵ = 0.02

  • D′: ϵ = 0.03

  • D_h: ϵ = 0.03.

• Minimum points: The minimum number of points required to form a core point in the cluster. This was set to 3, meaning each core point had to be associated with at least two other points.

The decoherence lifetime was defined as the time it took for the largest cluster to reduce to less than 5% of its original size, indicating significant dispersion of the meteoroid stream.

Fig. 1 Orbital distribution of 50 recovered meteorite falls (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.5. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

Fig. 2 Orbital distribution of 35 012 NEOs (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.05. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

Fig. 3 Orbital distribution of 310 US Government sensor impacts (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

3 Results

3.1 KDE fitting

In this study, we applied KDE to estimate the sporadic probability density functions (PDF) of the orbital distributions within our datasets. This approach is integral to assessing the statistical significance of orbital similarities, which will be explored in the next section. The KDE method was chosen because it effectively smooths the data to create a continuous representation of the underlying orbital distributions, possibly removing any small streams present in the data, while preserving the overarching features of each dataset. The PDFs (seen as red in the figures) estimated by the KDE method are shown for the 50 meteorite falls (Fig. 1), the 35012 NEOs (Fig. 2), 350 potential 50 g falls detected by GFO-FRIPON-EFN (Sect. 3.1), and the 310 USG Sensor detections (Fig. 3).

For the KDE fitting to be effective, it is crucial that the KDE accurately represents the general distribution of each dataset while smoothing out minor features such as small meteoroid streams or showers. This method is particularly suitable for datasets where showers constitute a minimal portion of the population. As demonstrated by Shober & Vaubaillon (2024), the KDE approach may not be applicable when showers make up a significant fraction of the dataset, as their presence can distort the overall distribution. To address this, fireball datasets with significant shower components need to have these showers removed prior to analysis. This can be achieved using a dissimilarity function that does not directly rely on the orbital elements, such as the Valsecchi et al. (1999) function.

In our datasets used here, the fireball data included only events with a high likelihood of dropping meteorites, effectively eliminating all shower contributions. The degree of smoothing in the KDE is controlled by the bandwidth parameter. For datasets with limited data, the bandwidth is set higher to prevent overfitting and ensure that any small streams or showers present are adequately smoothed out. This approach allows us to detect and account for any minor clusters within the dataset, ensuring that the KDE fitting accurately reflects the broader orbital distribution.

We utilized the KernelDensity class from scikit-learn to execute random sampling based on a Gaussian KDE (Pedregosa et al. 2011). The sampling process begins by randomly choosing base points from the dataset involved in the KDE fit, allowing each data point an equal probability of being selected. If sample weights are specified, these probabilities are accordingly adjusted. Each base point chosen, denoted as x_i, is then perturbed by adding Gaussian noise, thereby generating a new sample point x′_i calculated as x′_i = x_i + N(0, h²). Here, N(0, h²) represents Gaussian noise with zero mean and variance h², where h is the KDE bandwidth parameter. This parameter is crucial as it modulates the spread of the generated samples around their corresponding base points, thereby influencing the smoothness of the resultant density estimation. Each iteration of this procedure ensures that the sampled points closely mirror the density profile estimated by the KDE, using the Gaussian kernel’s properties to facilitate tasks such as Monte Carlo simulations and the generation of synthetic datasets. Detailed implementation specifics of scikit-learn’s KernelDensity class can be accessed through its GitHub repository⁹. This methodology of using KDE to estimate the random association rate between orbits was first employed by Shober & Vaubaillon (2024).

This random sampling from the estimated PDFs produced through KDE notably differs from previous studies that explored the chance similarity of orbits. Notably, the study of Jopek & Bronikowska (2017) primarily recommended the use of cumulative probability distribution (CPD) inversion to generate random samples by mapping uniform random variables to quantiles of the target distribution using the inverse of its cumulative distribution function (CDF). This method ensures that the generated samples accurately reflect the overall distribution’s quan-tiles. In contrast, KDE smooths observed data to estimate the PDF, providing a continuous approximation that highlights the density’s shape and local features while preserving parameter correlations. If some of the local features are real clusters and the CPD inversion includes them as part of the PDF to estimate random associations, the statistical significance of the cluster will be underestimated as the cluster was assumed to be a feature of the sporadic background. Meanwhile, the KDE’s ability to smooth out the possible presence of small clusters in the dataset when estimating the PDF removes this underestimation. With that in mind, choosing a bandwidth that is too large for the KDE and over-smoothing the PDF could equally result in an overes-timation of the significance of a cluster. Thus, regardless of the method chosen, one must be aware of these inherent setbacks and strengths.

In Pauls & Gladman (2005), they constructed two randomized orbital distributions by using the debiased NEO orbital distribution published by Bottke Jr et al. (2002) and from a dataset of 481 fireballs. They randomly selected (a, e, i) triplets from these distributions and used a pseudo-random number to determine if the selected triplet would impact Earth (only for NEO model). They then selected random longitudes of the ascending node uniformly from 0 to 360° and fixed the argument of perihelion such that either the ascending or descending node would be within the torus swept out by the Earth. This allowed them to construct databases of orbits with randomized angular elements but realistic semi-major axis, eccentricity, and inclination values. Albeit, this method also could suffer from the same slight underestimation of small clusters as discussed possible for the method of Jopek & Bronikowska (2017).

As seen in Fig. 1, due to the low statistics (only 50 orbits) of the meteorite fall with orbits dataset, the KDE’s smoothing can reduce the small dataset’s considerable variation and produce a more reasonable continuous distribution to draw from. The larger bandwidth of 0.5 used for the meteorite falls seems to produce an adequate distribution and is larger due to the small size of this specific dataset. The known NEO distribution comprises 35 012 telescopically observed objects (Fig. 2). Meanwhile, the impact database detected by USG Sensors has 310 events (Fig. 3) with velocity information available. Based on thousands of observations from the FRIPON, EFN, and GFO networks, 615 possible > 1 g meteorite-dropping fireballs were detected (Fig. 4), of which 315 fireballs are possible > 50 g meteorite-droppers (Fig. 5).

Fig. 4 Orbital distribution of 616 possible meteorite falls detected by GFO, FRIPON, or EFN sensors (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15 and ensuring any sample crossed the torus traced out by the Earth’s orbit. The meteorite falls were identified using α − β parameters as described in Sansom et al. (2019) using a minimum end mass of 1 g. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

Fig. 5 Orbital distribution of 350 probable 50 g meteorite falls detected by GFO, FRIPON, or EFN sensors (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15 and ensuring any sample crossed the torus traced out by the Earth’s orbit. The meteorite falls were identified using α - β parameters as described in Sansom et al. (2019) using a minimum end mass of 50 g. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

3.2 Statistical significance of orbital similarity

3.2.1 Fireball clustering

In Pauls & Gladman (2005), they found that despite previous claims, the probability that a chance association as similar as Příbram and Neuschwanstein (D′ ~ 0.009) was expected to occur ≅ 70% of time for a fireball dataset of 481 orbits. Even using another fireball model based upon a model NEA distribution Bottke Jr et al. (2002) still gave a chance occurrence of a D′ = 0.009 to be 10%, still not low enough to reject the null hypothesis. Furthermore, the other streams suggested by Halliday et al. (1990) were not identified as significantly lower than expected from simple random associations.

In this study, we have repeated this statistical analysis, using a similar KDE-based method to estimate the false positive rate as done in Shober & Vaubaillon (2024). Our model of the fireball population was fitted to 350 identified potential meteorite-dropping fireballs observed by the EFN, GFO, and FRIPON. This subset was selected through a combination of several metrics and methods, including the α − β methodology of Sansom et al. (2019), the fragmentation and ablation models used in Borovička et al. (2022b), shower identification and removal (Valsecchi et al. 1999), and sporadic cometary component removal (Tancredi 2014; Shober et al. 2021; Shober et al. 2024b). By randomly drawing 481 orbits from the PDF generated from this meteorite dropper distribution, we find that there is a ~3.1% chance of a random association. Despite this estimate being lower than that of Pauls & Gladman (2005), in the nearly 20 years since their study, there have been no more associations of that degree within the fall dataset. Thus, the likelihood of a chance association is higher than this, considering the many new and ongoing global fireball networks; the number of observed falls is very likely to be over a thousand or more. More importantly, we find that using our KDE false positive method, the meteorite-meteorite similarities using all three D-functions are within a 3-σ range expected for random associations (Fig. 6). The Příbram-Neuschwanstein pair only seems less likely for the D′ discriminant, while the other two metrics produce distributions that are very consistent with a spurious association.

The orbital distribution for the 50 recovered falls (Fig. 1), however, seems to significantly deviate from the distribution of predicted meteorite falls observed by fireball networks (Figs. 4 and 5). Thus, this distribution either (1) does not estimate well the population, or (2) searching and finders bias is significantly modifying the orbital distribution of recovered falls. When only using the 50 recovered fall data points to reconstruct the fall distribution, the statistical significance increases, at least when using a bandwidth of 0.5. However, using this KDE method on such a small dataset, while possible, is neither ideal nor trivial in selecting the bandwidth parameter well. It is much easier to accidentally over- or under-smooth the true distribution, leading to less precise results and possible false detections or false non-detections of streams. Thus, a larger dataset of 350 possible >50 g meteorite falls was used to estimate the fall population more precisely (Fig. 6), and the resulting correspondence between the observations and random association rate is considered more robust. Moreover, another critical factor is that there is a more general similarity amongst all possible pairs of the 50 recovered meteorites, as supported by the presence of a significantly shallower cumulative D-value slope, supporting the hypothesis that a search/recovery bias is affecting this recovered fall sample.

As explained in further detail in Section 2, our Monte Carlo KDE sampling method enabled us to fully explore the orbital similarity within the fireball, fall, and NEO datasets. The cumulative D-value distribution figures produced from this method demonstrate how all unique pair similarities within a single population or between two populations compare to those expected from random associations. The blue shaded areas in Figures 6–11 indicate the 3-σ regions given by our model consistent with the range of values produced by random associations. Meanwhile, the red lines indicate the actual observed cumulative D-value distribution.

First, we compared the similarity significance between the 50 meteorite falls and the known NEO population. Nearly every study of recovered meteorite falls with precisely measured orbits has tried to find NEOs that are in similar orbit space. Despite problems dealing with CRE ages and decoherence lifetimes, which we will address in Section 4, are any of these associations statistically significant?

These results are illustrated in Fig.7. The Monte Carlo simulation aimed to estimate the number of random associations expected for each of the three D-functions used here (D_sh, D′, D_h) using a dataset of 350 possible 50 g meteorite falls observed by the GFO, EFN, and FRIPON to estimate the orbital distribution. As depicted in the figure, the expected number of random associations closely matches the actual data, indicating that the orbital similarity between the meteorite databases (50 meteorite falls) and the NEOs is entirely consistent with random chance associations. This finding contradicts several studies that have used low D values with NEOs to suggest associations between meteorites and near-Earth objects. Similar results were also found when the meteorite fall PDF was estimated using the orbits from the 50 recovered meteorites population (Fig. 1), and the 616 possible 1 g falls (Fig. 4). No matter the fireball subset used, we found no statistically significant clustering between the fall and NEO populations. Our results demonstrate that the associations are consistent with random chance throughout the D-value ranges. One interesting aspect, however, is the noticeable overestimation of the KDE at larger D values. We also tested the random association likelihood drawing from a KDE fit to a very smoothed (bandwidth=0.5; Fig. 1) distribution of the 50 recovered meteorite falls, and this also produced cumulative D-value distribution well within the 3σ region of random associations.

D values have also been used to identify the parent or associated asteroids from the NEO population with USG sensor-detected impacts. Thus, we also compared this population with the NEO population to estimate the number of random associations as a function of the D value. Fig. 8 shows that the data is consistent with random associations for all three metrics. This analysis indicates no evidence of statistically significant similarity between the CNEOS database and the NEO population.

The CNEOS bolide database, which uses USG sensor observations, is widely known now to have substantial and seemingly unpredictable orbital precision, and no formal uncertainties are associated (Devillepoix et al. 2019). On the other hand, for fireballs observed by dedicated photographic networks, where the orbits and their corresponding uncertainties are understood, can we find statistically significant clustering? Using a lower final mass limit of ~1 g to augment the size of the dataset, we collected 616 fireballs from GFO-FRIPON-EFN observations (Fig 4). Despite the increased precision of this data, as seen in Fig. 9, the number of associations is also completely in correspondence with that expected from random associations. For the moment, there seems to be no evidence based on orbital similarity that streams exist in the centimeter- to meter-sized populations that are impacting the Earth’s atmosphere.

In the analysis of fireball and meteorite recoveries, there is often a logical assumption that low D values are less significant for bodies with inclinations near the ecliptic due to the high density of such orbits. In contrast, similar D values for higher inclination orbits are considered more meaningful. In this study, we have quantified how the statistical significance of these D values varies as a function of the orbit as well. Thus far, we have only examined whether the total number of orbital pairs is below a certain threshold. However, we must also consider the probability of getting an orbital pair below some D-threshold within the context of the orbital space of the pairs to avoid the missed identification of small clusters in low-probability orbits (e.g., high inclinations). It has also been shown recently how theoretical NEO families disperse more slowly when they are formed at higher inclinations and further from the nearest planet (Humpage & Christou 2024). Thus, any clusters on these orbits should persist longer, and the examination of the statistical significance as a function of the orbit cannot be overlooked.

Nevertheless, as seen in Fig. 10, we find that even when the statistical significance is tested as a function of the orbit, we still find no pairs that exceed a 3σ limit. To construct Fig. 10, we employed the same methodology as used in creating Fig. 9 using the data of the 616 GFO-FRIPON-EFN fireballs. However, we also divided the orbital space into bins, and for each bin, we tracked the number of pairs with D_H < 0.1. The heat map coloration in Fig. 10 then shows the mean number of such pairs over all the Monte Carlo trials, plus three standard deviations. This approach allowed us to visualize where, within the orbital space, the statistical significance might deviate from random expectations. We also tested this mapping for D_h < 0.05 and D_h < 0.1 and logically encountered the same results; in other words, no statistically significant groupings.

Our findings reveal that none of the fireball pairs within this dataset exceed the 3-sigma limit, indicating no statistically significant clustering, even when accounting for the specific orbital contexts. We find only one pair (red point) amongst the 616 possible 1 g falls that satisfies the D_H < 0.1 used to make the heat map in Fig. 10, but it does not surpass a 3σ significance threshold. Additionally, we re-considered the pair of Příbram-Neuschwanstein (black square) in Fig. 10. However, this too also failed to meet a 3σ significance threshold when considering a dataset of 481 orbits, as done in Pauls & Gladman (2005). A similar analysis of the CNEOS dataset identified one impact event above the 3-sigma threshold at a 10° inclination. However, this result becomes inconsequential when considering the inherent unconstrained orbital uncertainties in the CNEOS dataset. Another similar claim was made for the Chelyabinsk impactor and NEA 1999 NC43; however, based on compositional and spectral information, this link was considered unlikely (Reddy et al. 2015). Beyond this non-correspondence, the statistical significance test used previously has some pretty large implicit assumptions that we would like to discuss. The statistical significance was determined by calculating the probability of drawing a randomly selected asteroid that is as similar to the measured Chelyabinsk orbit, only comparing to the asteroids with larger diameters than 1999 NC43 (Borovička et al. 2013). To robustly reject this null hypothesis, one must assume that Chelyabinsk is not unique and is another random draw from the same source population of all other observed falls. This is the same problem as the “birthday paradox” discussed in Pauls & Gladman (2005). This veridical paradox states that within a group of 23 people or more, there is a greater than 50% chance of two people sharing a birthday. This is not a true paradox, but rather just an unintuitive result, as most people approach this problem by trying to calculate the probability that they will personally share a birthday with someone in the group. This is the same logic used when only testing the probability of similarity with a specific orbit instead of the probability of a random association between any of the combinations of random draws of the concerned population(s). Additionally, the limitation of the asteroid population to only bodies greater than the size of 1999 NC43 significantly increases the significance of the association without adequate justification of the orbital uniqueness of this larger asteroid subset. In this study, we have redone this significance test using our KDE-based methodology. From this we find that when using similar assumptions (Borovička et al. 2013; Reddy et al. 2015), similar significance values are found. However, considering that there were 20 recovered meteorites at the time of the Chelyabinsk impact and assuming this sample is being drawn from the same source population, the probability of a random association is still more than 3σ. However, if the size of the asteroid population includes smaller bodies, the orbital similarity is of 1999 NC43, and the Chelyabinsk impactor becomes insignificant. Thus, when calculating the statistical significance of two bodies, it is important to consider the assumptions carefully.

Thus, by quantifying the relationship between D-value significance and orbital parameters, we find no robust evidence of statistically significant clustering. This further supports the conclusion that observed similarities are consistent with random associations rather than true physical streams in impact datasets.

Fig. 6 Number of meteorite-meteorite pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and DH. The blue circles represent the mean number of pairs from 500 iterations of KDE sampling using a bandwidth of 0.1 based on 350 possible/likely 50 g meteorite falls observed by the GFO, FRIPON, and EFN networks. The blue-shaded area is the 3 σ range possible from random associations. The actual number of low-D value pairs is consistent with the predicted number based on random associations.

Fig. 7 Number of meteorite-NEO pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H , D′, and DH. The blue circles represent the mean number of pairs identified between 50 random samples of “meteorite falls” with the observed 35012 NEOs from 500 iterations ofKDE sampling using a bandwidth of0.15 forthe350possible50gmeteoritefalls (Fig. 5). Theblue-shaded area is within three standard deviations for random samples of 50 (the number of meteorites with orbits). The actual number of low-D value pairs correctly follows the predicted number based on random associations, indicating no evidence supporting the presence of streams in the recovered meteorite fall dataset. Two asteroids were removed as they correspond to two recovered meteorite falls, Motopi Pan and Almahata Sitta, observed telescopically before impact.

Fig. 8 Number of CNEOS USG Sensors impact – NEO pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and DH. The blue circles represent the mean number of pairs identified between 310 random samples drawn from the USG impactors PDF (Fig. 3) with the observed 35 012 NEOs from 500 iterations of KDE sampling. The blue-shaded indicates a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations using the KDE distribution shown in Figs. 2, indicating no substantial evidence supporting streams in the dataset.

Fig. 9 Number of unique pairs, ND, within the 616 possible 1 g final mass fireballs observed by the GFO-FRIPON-EFN with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and Dh. The blue circles represent the mean number of pairs identified between two 616 random samples drawn from the possible 1 g fall PDF (Fig. 4) from 500 iterations of KDE sampling. The blue-shaded indicates a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations, indicating no substantial evidence supporting streams in the fireball dataset.

Fig. 10 Maximum false positive heat map for bins of semi-major axis (au) by inclination (deg). The heat map coloration represents the mean number of pairs plus 3σ that have DH < 0.1. The sporadic orbital distribution was estimated using the 616 possible 1 g fall observations PDF (Fig. 4). The means and standard deviations were calculated through 500 iterations of KDE sampling. Subplot a examines the maximum number of expected random associations within a population of 616 fireballs, only one (red) in the observations has a measured DH < 0.1, but it is not statistically significant. Subplot b re-examines the significance of the Príbram-Neuschwanstein pair (black square); however, the pair is still not statistically significant within a fireball population of 481 orbits.

Fig. 11 Number of NEO – NEO pairs with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for Dsh, D′, and Dh. The blue circles represent the mean number of pairs from 500 iterations of KDE sampling using a bandwidth of 0.15, with a shaded blue region indicating a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations for large values; however, at very low D values, there appear to be tens of objects with statistically significant similarities. The observational biases of the dataset do not predict the elbow in the cumulative D-value plots.

3.2.2 NEO clustering

Several investigations have also explored the presence of NEO families (Fu et al. 2005; Schunová et al. 2012; Jopek 2020). Using our KDE-based methodology, we have re-examined the statistical significance of orbital similarity among all unique pairs within the 35 012 NEOs from the NASA HORIZONS database. Our analysis revealed a remarkable consistency in the distribution of large D values across all three orbital similarity discriminants (D_{s h}, D_h, and D′). This is clearly illustrated in Figure 11, where the observed cumulative D-value distribution shows a distinct deviation from the random association model, particularly around D_h = 0.03. This deviation, or “kink,” which is present across all similarity distributions, indicates the presence of statistically significant groupings within the NEO population – a feature that contrasts sharply with the expected results from random associations. This kink is also not completely explained due to the observational bias of telescopic surveys, as that would have been captured in the blue 3-σ region in Fig. 11. The statistical significance of these groupings begins to emerge at D values around 10⁻², consistent across the three discriminants, although with some variation. Previous studies, such as Schunová et al. (2012), found no consistent evidence for NEO families, while more recent work by Jopek (2020) identified 15 potential groupings with over 900 objects. Our results corroborate these more recent findings that statistically significant clustering is occurring, albeit we find that the clusters identified in Jopek (2020) are much more inclusive than we find here.

We used the DBSCAN clustering algorithm employing D_H, an epsilon value of 0.03, and a minimum of two connections for a pair to be defined as a core point. This revealed clusters that tend to align near Earth-crossing orbits, as shown in Fig. 12. The epsilon value was chosen precisely where we see a “kink” in the cumulative D-value distribution (Fig. 13), indicating that a different production mechanism becomes dominant. This is just slightly larger than where we find the statistical significance is beyond 3-σ for D_h (~0.01), but significantly less than the minimum D_H values used in Jopek (2020). Jopek (2020) also uses a single-linkage clustering approach, which tends to produce elongated or “chain-like” clusters, especially in the presence of noise. The difference in clustering algorithm and larger D-value limits are likely the primary reasons why Jopek (2020) identified hundreds of clustered NEOs compared to the tens that we identify here. That being said, we have not estimated the probability of having a group of pairings like done in Jopek (2020); thus, we are only plotting a minimum number of clusters where each linkage is statistically significant.

In Fig. 12, we identified 12 clusters of at least 5 members satisfying our DBSCAN identification scheme. We found 66 clustered NEAs distributed between the 11 clusters, 28 of these objects being previously identified by Jopek (2020). These new identifications are mostly due to the significant increases in the NEO dataset size over the previous two decades. The NEO dataset used by (Jopek 2020) consisted of 20032 NEAs, Schunová et al. (2012) analyzed 7563 NEOs, and the study by Fu et al. (2005) only had 3319 NEOs to compare. Thus, the significantly larger dataset of 35 012 NEOs considered here is more amenable to discovering small ephemeral NEO families. However, this increase in NEO discoveries does come with its caveats, as surveys do not observe all objects equally.

In Fig. 12, the heatmap shows the observational bias by comparing the synthetic debiased NEO dataset of Granvik et al. (2018) with the observed NEO dataset (35 012 objects). The red regions represent an overabundance of NEOs due to a discovery bias, while blue areas indicate under-representation in observations. Crucially, while we have demonstrated through our KDE Monte Carlo method that there is a statistically significant number of pairs on similar orbits as shown in Fig. 11, all the clusters identified are in regions where an overabundance of NEOs is observed due to this discovery bias, as highlighted in Fig. 12. This does not necessarily imply that these clusters are entirely a result of the discovery bias, as our KDE Monte Carlo method should account for such biases, ensuring that the observed clustering is indeed statistically significant. However, the number of clustered objects identified by Jopek (2020) is could be an over-estimation based on our results, and this is possibly influenced by this observational bias. Notably, as seen in Fig. 12, there are regions, particularly on orbits with semi-major axes between 0.5 and 1.5 au and higher eccentricities, where one would expect to see more clusters if they were solely due to observational bias. The absence of clusters in these regions supports our finding that while observational bias likely slightly contributes to the clustering identified, it cannot account for all of it. The clustering that we have identified also successfully includes the fragments of the well-observed breakup of comet 73P/Schwassmann-Wachmann (yellow group; Fig. 12) as a statistically significant cluster, distinct from the background, without prior intent for identification, giving us confidence that these clusters were identified using adequate constraints.

This clustering observation is also consistent with recent studies suggesting that tidal disruptions are responsible for creating such groupings (Granvik & Walsh 2024). Our analysis identified approximately 10^1–2 such objects according to Fig. 13, supporting the hypothesis that tidal disruptions contribute to the observed clustering. These results interestingly contrast with orbital meteorites, fall data, and the CNEOS impact data set, for which no statistically significant pairings exist. This disparity between the NEO groupings and the meteorite and impact data remains an intriguing point for further investigation, particularly in understanding the differing influences of tidal disruptions and other dynamic processes on these populations. Other mechanisms, such as rotational disruptions, meteoroid impacts, etc., are other possible debris generation mechanisms worthy of exploring, particularly for smaller debris (Jewitt 2012).

Given the clustering observed in the NEO dataset, which accounts for up to a few percent of the population (Jopek 2020), a similar proportion of clustering would be expected in the fireball and impact datasets if these populations are comparable. To explore this hypothesis, we also integrated the fireball data with the NEO dataset and reapplied the DBSCAN clustering algorithm utilized in generating our NEO clusters Fig. 12. Surprisingly, this combined analysis revealed no changes in the identified clusters; the same clusters persisted in the NEO dataset without any inclusion of fireball data. This unexpected result suggests two possibilities: either the inherent uncertainties in fireball data are effectively smoothing out potential small associations, or the proportion of clusters within the impact dataset is actually lower than what we observe within the NEO population.

This decrease in the proportion of detectable clusters could also be influenced by the intrinsic biases of the two datasets. The NEO dataset and the impact dataset are both biased, but in opposite ways. The NEO dataset is skewed towards detecting larger and brighter objects on orbits near the Earth (Fig. 12), which tend to have lower encounter velocities and more Earth-like orbits. On the other hand, the fireball dataset is biased towards higher velocity impacts, particularly those from smaller objects, because higher velocity impacts tend to be brighter and thus more likely to be observed. Consequently, this dataset is more likely to capture events with larger semi-major axes, higher eccentricities, and greater inclinations, the opposite of the bias seen in the NEO dataset. Therefore, these conflicting biases could also contribute to the lack of observed associations between the fireball data and the NEO clusters.

To extend this, we have also compared the observed cumulative D-distributions (i.e., the distribution of all distance values between all objects in the population(s)). This was visualized in a log-log plot in Fig. 13, where the y-axis represents the number of objects with a D value less than that on the x-axis. We observed that datasets showing no evidence of statistically significant associations exhibited a linear distribution in log-log space for this cumulative plot. We calculated the orbital similarity of all unique pair combinations within the USG sensor data, the 50 recovered meteorite falls, 616 fireballs observed by GFO-FRIPON-EFN, and the sporadic subsets of the FRIPON and EFN datasets after removing known meteor showers using the Valsecchi D-function (Valsecchi et al. 1999). This function was only used to remove the shower component as it does not rely on orbital similarity directly, avoiding the self-fulfilling bias of the cumulative D-value slope analysis. The linearity in log-log space suggests a lack of significant streams, as there is no excess of pairs with small D values to indicate clustering beyond what is expected by random chance. This linear behavior can be explained by the random distribution of orbital elements. In contrast, datasets with significant streams show a deviation from this linear trend, with a change in slope at smaller D values, indicating an excess of pairs with high orbital similarity. This is demonstrated clearly through the contrast between the FRIPON and EFN complete datasets and the sporadic sub-component.

Our analysis indicates that the populations without evidence of streams, including the USG sensors and the 50 recovered meteorite falls, follow this linear trend for a cumulative D-value distribution. This provides further evidence that no significant meteoroid streams are present in these datasets. The lack of deviation from linearity in the log-log cumulative plots confirms that the observed orbital similarities are consistent with random associations, reinforcing the conclusions drawn from our KDE-based Monte Carlo simulations. The cumulative D-value distributions for the sporadic populations exhibit very similar slopes, suggesting a general degree of similarity between their distributions and, thus, self-similarity. Notably, the 50 meteorite falls deviate from this pattern by displaying a slightly lower slope while maintaining linearity. This lower slope does not resemble the FRIPON and EFN datasets with showers included. Additionally, it contrasts with the observed fall database of 616 potential meteorite-dropping fireballs observed by the GFO, EFN, and FRIPON. This difference likely indicates a discovery or recovery bias influencing our recovered meteorite database. This is not entirely surprising, given the challenges and limitations of meteorite recovery.

Notably, the NEO dataset, as previously discussed, does display a significant slope-change or “kink” in the cumulative D-value distribution around 10⁻². This clear change in slope supports the hypothesis that tidal disruption or another mechanism in near-Earth space generates short-lived small NEO families (Schunová et al. 2014; Granvik & Walsh 2024). However, the absence of a significant change in slope in the cumulative D-value plots across all other analyzed datasets corroborates our findings that there is no substantial evidence for meteoroid streams among meteorite falls or CNEOS impacts.

Fig. 12 NEO orbital distribution with ephemeris taken from NASA Horizons with clusters identified, objects not belonging to a cluster can be seen as small gray points. The background heatmap, however, indicates the observational bias of the NEO dataset. The red regions denote orbits for which more objects exist in the NEO database relative to the debiased dataset of Granvik et al. (2018), and blue is where there exists a deficit of NEO detections. Clusters identified by Jopek (2020) can be seen as large dark gray points. The clusters identified within this study are colored, each with at least 5 members. They were identified using a DBSCAN algorithm where core points were defined as having at least two associations with a minimum ϵ corresponding to DH value of 0.03. This D value was chosen as it corresponds well to where the cumulative D-value distribution displays a “kink” corresponding to an excess of similarity. The points that were also identified by Jopek (2020) as being in a cluster have a black border. The new large yellow cluster identified is not asteroidal but corresponds to all the fragments of 73P/Schwassmann-Wachmann, which was a Jupiter-family comet observed to have undergone significant fragmentation between 1995 and 2006 (Reach et al. 2009).

Fig. 13 Cumulative D-value distributions for every possible pair within 310 USG Sensors orbits 3290 FRIPON fireballs, 824 EFN fireballs, 1913 FRIPON sporadic subset fireballs, 457 EFN sporadic subset fireballs, 50 recoveredmeteoriteorbits, 616potential1 gmeteoritedroppers observed by GFO-FRIPON-EFN, and 35 012 NEOs taken from the NASA Horizons ephemeris service.

3.3 Lyapunov characteristic lifetimes

Building upon the analysis of orbital similarities and the identification of statistically significant clusters within the NEO population, we now turn our attention to the underlying dynamical behavior of these objects. The chaotic nature of near-Earth space, influenced heavily by gravitational interactions, suggests that the orbital paths of these objects are not only influenced by their current configuration but are also highly sensitive to initial conditions. To explore this, we investigated the Lyapunov characteristic lifetimes of NEO orbits, which quantitatively measure their dynamic stability. Specifically, we examine how these lifetimes vary with orbital parameters such as semi-major axis, eccentricity, and inclination. We focus on their implications for the temporal evolution of meteoroid streams and the practical limits of using orbital similarity metrics like the D values over extended periods.

Figure 14 presents the Lyapunov characteristic lifetimes as functions of the semi-major axis and eccentricity for three distinct inclinations: 0°, 10°, and 20°. These maps were produced using the Rebound N-body simulation software. The model of the solar system included the Sun and the eight major planets, utilizing initial conditions sourced from JPL Horizons. The ‘whfast’ integrator, known for its efficiency with near-Keplerian orbits, was employed with a time step set at 0.01 years (Rein & Tamayo 2015).

A dense grid of points represents initial orbital parameters for hypothetical small bodies for each map with linear interpolation in lower statistic regions. Specifically, 289 640 particles were integrated at 0° inclination, 287 352 at 10° inclination, and 353 470 particles at 20° inclination, covering a semi-major axis range from 0.5 to 3.5 au and eccentricities from 0 to 0.95.

The color scale in each map indicates the Lyapunov lifetime, which quantifies the rate at which two infinitesimally close orbits diverge. This highlights the system’s sensitivity to initial conditions. The simulations were run over a period of 20 000 years, and the Lyapunov time was calculated as the inverse of the Lyapunov exponent.

The Lyapunov lifetime maps exhibit significant variation across different orbital elements and inclinations. At a low inclination of 0°, the maps display considerable variability in lifetimes, with shorter times noted in regions of higher eccentricity as they cross the orbits of more planets. Notable, too, are mean-motion resonances such as the 3:1 and 2:1 with Jupiter, which appear as zones of enhanced stability on these timescales (longer Lyapunov lifetimes). This is linked to close encounters: meteoroids inside them are unable to encounter the planet they are in resonance with. With an increase to a moderate inclination of 10°, the chaotic zones tend to be further concentrated towards the aphelion and perihelion lines of the planets as these bodies tend to have close encounters. At the highest studied inclination of 20°, chaotic zones are more widespread, especially for orbits intersecting multiple planetary paths. The protective effect of mean-motion resonances lessens, and overall stability across the semi-major axis range decreases, suggesting that higher inclinations generally lead to more chaotic orbits in near-Earth space (except for more circular orbits). This protection of the mean-motion resonance is particularly not visible for large eccentricities when the inclination reaches 20° . A previous work (Courtot et al. 2024) has shown that a chaos indicator in that region reveals regular dynamics. This could be explained in two ways: this area could be qualified as “stable chaos,” or the chaotic aspect could also appear only on longer timescales, such as the 20 000 years used here.

Comparing the three maps with each other, the evolution of the Lyapunov lifetime as the inclination rises seems to be strongly linked with close encounters. As the meteoroids are further away from the ecliptic, the zones of shortest Lyapunov time change shape and follow less clearly the zones of close encounters delimited by the black lines. It is well-known that the geometry of close encounters, and thus its effects on meteoroid streams, depends on the direction of the incoming meteoroid with respect to the planet’s orbit (Carusi et al. 1990). This explains how the shape of areas of shorter Lyapunov lifetime can be distorted for higher inclinations. We can also note that for the inclination at 10° , therefore close to Mercury’s orbital inclination at 7.004°, Lyapunov lifetimes are especially short when close encounters with Mercury are possible, meaning Mercury would be the driving force behind the chaos in this configuration.

For the inclination at 0° , a large portion of the map includes Lyapunov lifetimes shorter than 100 years, while this value can go down to 60 years for a large part of the map for the inclination at 20°. The map at10° shows less areas with shortLyapunov lifetime, except for encounters with Mercury, as discussed above. Therefore in those regions, meteoroids with initially similar orbits might diverge quickly, rendering the use of D-functions quite complex. A short Lyapunov lifetime reveals areas where meteoroids on initially similar orbits might reach widely different orbits. Even if D values are able to characterize whether the orbits are similar today, it does not predict how they will evolve: an existing stream today might disappear quite fast. Therefore, a short Lyapunov lifetime in an area only shows where any D-function results must be taken with a grain of salt.

Finally, it should be noted that we are especially interested in the area where meteoroids can encounter the Earth, since D-fucntions are mostly used to find meteor showers. Taking this condition into account, the Lyapunov lifetime is usually very short (less than 200 years for i = 0°, less than 500 years for i = 10° and less than 100 years for i = 20°), and it gets much longer on only very small areas. This means that for the conditions we are especially interested in, results from D-functions should be regarded with caution.

The findings from the Lyapunov lifetime maps have significant implications for the interpretation of statistical clusters in near-Earth space, as discussed in the previous section. Despite the presence of tidal disruptions in near-Earth space, the orbital paths that result from these disruptions also tend to be the most chaotic regions. The clusters identified in Fig. 12, correspond to Lyapunov characteristic lifetimes of only tens to hundreds of years – as the perihelia values tend towards 1 au. Any clustering of objects resulting from recent tidal disruption events would have a very short dynamical memory. Consequently, while we can statistically analyze orbital similarities among near-Earth objects, precisely identifying which objects belong to a particular cluster becomes problematic due to the rapid divergence of orbits within these chaotic regions.

This challenge is further underscored by the contrast between the frequency of tidal disruption events and the typical Lyapunov lifetimes. For instance, the estimated frequency of tidal disruptions with the Earth that would result in a detectable cluster in the NEO population has been estimated previously to be anywhere between a few to several thousand years (Richardson et al. 1998; Tóth et al. 2011; Schunová et al. 2014). Thus, while tidal disruption could generate clusters, the chaotic nature of their orbits means that meaningful connections between the bodies in these clusters would disperse very quickly. In conjunction with higher intrinsic uncertainties of orbits extracted from short observation arcs as a body impacts the atmosphere, this chaotic dispersion may also help explain why we do not observe statistically significant clustering in fireball datasets. With increased fireball data, clusters may become statistically significant, but perhaps this will require thousands to tens of thousands of very precise fall observations. Currently, the lack of statistically significant clustering in fireball data suggests that any observed similarities are more likely to be spurious associations rather than evidence of real clustering. Any current claims linking specific meteorites to near-Earth objects should be viewed skeptically. The statistical analysis presented here indicates that such associations are not statistically sound and are more likely to result from random chance rather than genuine physical connections.

Fig. 14 Lyapunov characteristic lifetime as a function of the semi-major axis and eccentricity with an inclination of 0° (a), 10° (b), 20° (c). The black lines indicate where the aphelion and perihelion distances equal the nominal semi-major axis value of Mercury, Venus, Earth, Mars, or Jupiter.

3.4 Meteoroid stream decoherence lifetimes

While the Lyapunov characteristic lifetime offers valuable insights into the chaotic nature of individual meteoroid orbits, it does not fully encapsulate the persistence of coherent structures within a meteoroid stream. The Lyapunov lifetime typically reflects the timescale over which small perturbations lead to exponential divergence in the orbits of individual particles, usually on the order of a few hundred years in near-Earth space. This measure primarily addresses the evolution of one particle compared to others with similar initial conditions rather than the longevity of the collective behavior of a meteoroid stream. In contrast, the decoherence lifetime of a meteoroid stream focuses on the collective dynamics of the stream’s constituent particles, quantifying the duration over which the stream remains identifiable as a distinct, statistically significant entity, despite the underlying chaotic dynamics. This lifetime often spans tens of thousands of years, significantly longer than the Lyapunov lifetime, because it considers the overall correlation in the orbits of the group of meteoroids, which can remain similar enough to maintain the stream’s identity over extended periods.

A potentially helpful analogy to illustrate this difference is to think of a cloud of smoke rising from a chimney. As the smoke particles ascend, small air currents cause them to move in different directions, leading to the rapid dispersion of individual particles. The point at which the paths of individual smoke particles start diverging from each other is analogous to the Lyapunov lifetime, where minor differences in initial conditions lead to rapid divergence in the orbits of meteoroids. However, despite the individual particles dispersing, the cloud itself remains visible as a cohesive entity for a considerable time as it drifts away from the chimney. The overall shape and presence of the smoke cloud linger, even as individual particles follow chaotic paths. This extended visibility and coherence of the cloud represent the decoherence lifetime, during which the smoke continues to be recognized as a collective stream despite the chaos experienced by individual particles. This analogy helps to clarify how a stream can remain a discernible group long after its constituents have started to exponentially diverge.

The study by Pauls & Gladman (2005) estimated the decoherence lifetimes for the three hypothetical meteoroid streams – Příbram, Innisfree, and Peekskill – also by generating test particles with nearly identical Keplerian orbits. These particles were integrated 500 kyrs considering the gravitational influences of all nine planets. However, since they only considered three streams, they did not need to place a strict condition on when decoherence officially occurred. They defined the stream as decoherent when the orbits of the particles had spread out sufficiently that any tight clustering in their orbits, which would have indicated a coherent stream, had disappeared. This loss of coherence was monitored by observing the dispersion of the orbits over time, specifically by examining how the encounter orbital longitudes of particles with Earth spread out.

In our study, we sought to define the decoherence lifetimes of meteoroid streams in a manner that would allow us to systematically compare the persistence of these streams across different regions of orbital space. This was important as a more qualitative definition of decoherence was impractical, given that we had generated 300 fictitious streams in total to achieve this goal. To do this, we considered the stream to have become “decoherent” when the largest remaining cluster within the stream only constituted 5% or less of the original size of the stream. This threshold was selected to ensure that we were capturing the significant dispersion of the stream beyond mere perturbations, reflecting a point where the stream’s identity has effectively been lost. We once again used the DBSCAN clustering algorithm, as it allowed us to identify and track clusters over time.

The choice of ϵ values in our DBSCAN clustering analysis is critically important, as it directly influences how we define and detect the coherence of meteoroid streams. The ϵ value determines the maximum distance between two points to be considered part of the same cluster, effectively setting the sensitivity threshold for identifying coherent structures within the data. The key challenge here is that the appropriate ϵ value is not constant; it must be carefully adjusted based on the size of the dataset and the proportion of the dataset that is made up of meteoroid streams. When streams constitute a small fraction of the overall data, smaller ϵ values are required to avoid mistakenly merging distinct groups or failing to detect smaller streams. Conversely, larger ϵ values can be used in datasets where streams form a more significant portion to capture the broader structures and avoid fragmenting real streams into smaller, artificially distinct clusters.

This dependency is illustrated in Fig. 13, where we observe distinct changes in the slope of the cumulative D-value distributions for the FRIPON and EFN databases compared to the NEO dataset. The D values where these changes occur are approximately an order of magnitude larger for the FRIPON and EFN databases than for the NEO dataset. This difference is primarily due to the higher proportion of cometary meteoroids in the FRIPON and EFN data, which make up a significant fraction of the observed orbits. Specifically, as shown in Shober & Vaubaillon (2024), meteor showers constitute roughly 18–25% of the EFN dataset. In contrast, streams in the NEO dataset represent only a small fraction, possibly just a few percentage points at most (Jopek 2020). This disparity necessitates a larger ϵ for the FRIPON and EFN databases to cluster the significant stream components appropriately. In contrast, with its more dispersed and sparse stream population, the NEO dataset requires a much smaller ϵ to identify the subtle clustering signals.

Our specific choice of ϵ values was thus strategically tied to the underlying processes we aimed to study, given the specific populations at hand. The ϵ value for the DBSCAN algorithm was chosen based on the “kink” observed in the cumulative distribution of the D values for NEOs (Fig. 13). We interpret this unique slope change as a signature of tidal disruptions -a mechanism distinct from the cometary outgassing that dominates the FRIPON and EFN datasets. By aligning our ϵ selection with this kink, we targeted the detection of streams generated through tidal disruptions, allowing us to map the decoherence lifetimes of these streams within near-Earth space. This approach ensures that the resulting decoherence maps are representative of the decoherence lifetimes of streams generated from tidal descriptions or any other mechanism that is causing this excess of similarity amongst NEOs. We interpret this “kink” as being likely linked to tidal disruptions, as this phenomenon has been supported by recent studies (Schunová et al. 2014; Granvik & Walsh 2024). These are the same ϵ values as those chosen in the cluster identifications of Fig. 12.

The resulting decoherence maps, presented in Figs. 15, 16, and 17, illustrate how these lifetimes vary as a function of orbital space and the similarity criterion used. At zero-degree inclination, we observe that Earth-crossing orbits generally have shorter decoherence lifetimes, on the order of tens to hundreds of thousands of years. As the inclination increases to 10 and 20 degrees, these lifetimes extend but still remain within this same general range. These results are consistent with those found by Pauls & Gladman (2005), indicating shorter lifetimes for streams on Earth-crossing orbits. Furthermore, the influence of Jupiter is apparent, with streams at higher inclinations or those closer to Jupiter’s orbit experiencing more rapid decoherence due to the planet’s gravitational perturbations. This indicates that while tidal disruptions may create clusters, the chaotic nature of these orbits results in a rapid loss of coherence, emphasizing the ephemeral nature of such streams in near-Earth space.

There also appears to be some noticeable variation between the orbital similarity discriminants used in this study. While the overall trends in decoherence lifetimes as a function of semi-major axis, eccentricity, and inclination remain consistent across the different discriminants, the absolute values of the lifetimes do show some variation. For instance, the D_H and D_{S H} discriminants exhibit nearly identical distributions, with DS H showing slightly, though very minutely, longer decoherence lifetimes. In contrast, the D′ discriminant consistently shows longer deco-herence lifetimes, a difference that is likely attributable to its stronger dependence on eccentricity (Jopek 1993). This effect is evident in Figs. 15, 16, and 17, where the decoherence lifetimes associated with D′ are slightly shifted towards larger values compared to those of D_H and D_{S H} regardless of the inclination. D_H and D_{S H}, which are more dependent on perihelion distance than D′, show more rapid changes in orbital similarity due to these encounters, leading to shorter decoherence lifetimes. However, despite these differences, the overall trends remain consistent, and the decoherence lifetimes for Earth-crossing meteoroid streams across all three discriminants are generally within the 10–100 kyr range with evolved orbits (semi-major axis no longer in range of the main-belt) capable of longer decoherence lifetimes.

Fig. 15 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 0° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. The ϵ value for the DBSCAN algorithm was chosen to where the “kink” is in the cumulative D-value distributions for NEOs (Fig. 13). Also, core points have at least 2 connections. The decoherence lifetime was thus defined as when the fictitious meteoroid stream lost 95% or more of the original stream mass.

4 Discussion

Using a kernel density estimation (KDE) method modified from Shober & Vaubaillon (2024), we evaluated the statistical significance of orbital pairings within several meteoroid and near-Earth object (NEO) datasets. This method allowed us to assess over 3 × 10¹¹ possible pairs amongst terrestrial impact datasets and NEOs, supplemented by Monte Carlo simulations to rigorously quantify the significance of these pairings, as outlined in Section 2.

Our analysis revealed that within the meteorite fall database, comprising 50 meteorites with well-documented orbits, there is no conclusive evidence of statistically significant pairings. Furthermore, when compared to the NEO clusters, only 7 fireball events of the 616 possible 1 g falls observed by the GFO, FRIPON, and the EFN were identified to satisfy the DBSCAN clustering σ and N_min constraints as used for Fig. 12. None of the 50 meteorite falls or the 310 CNEOS impacts were identified to be clustered. The cumulative distribution of D values (D_{S H}, D′, and D_H) across these datasets was linear, which aligns with random chance associations rather than genuine clustering. Also, more importantly, the cumulative D-value distributions corresponded well with the number expected due to chance association. Even when considering clustering with the NEO population, only ~1.2% are identified to be similar to identified NEO clusters with 5 or more members. Despite numerous studies attempting to link specific meteorite falls with NEOs via orbital similarity, our findings suggest that all these pairings lack statistical significance. The linearity observed in the cumulative distribution of D values strongly indicates that any observed orbital similarities between meteorites and NEOs are coincidental and not indicative of actual streams.

In contrast, our analysis confirms the presence of statistically significant clusters within the NEO population, corroborating the findings of Jopek (2020). However, our results suggest a more conservative estimate of the number of such clusters, identifying 12 statistically significant clusters compared to the 15 reported by Jopek (2020). Despite this difference in numbers, the locations of these clusters – primarily around perihelion distances near Earth’s orbit and in low-inclination orbits – are consistent with previous studies by Schunová et al. (2014) and the recent findings by Granvik & Walsh (2024) for evidence of tidal disruption in near-Earth space. The clustering of NEOs in these regions aligns with the hypothesis that tidal disruptions of NEOs, particularly during close encounters with Earth and Venus, produce a statistically significant excess of objects in these orbits.

It is essential to note that this does not necessarily mean that no meteorites are being produced from this tidal disruption mechanism. Quite the contrary, in a very recent work, Shober et al. (2024a) argues that the release of CI and CM chondrites from an immediate precursor body already in near-Earth space would much better explain the CRE ages, petrographic features, mixture of irradiated grains, and isotopic geochemistry in the samples. There may be many reasons why no meteorite with orbits and fireball observations are above the detectable limit. Firstly, the tidal disruption mechanism is not so dominant in near-Earth space that it produces a significant proportion of the debris. According to Jopek (2020), only 4.7% of the NEOs are identified in clusters, whereas we find this value to be lower, possibly below 1%. In any case, these minor streams require significant amounts of meteorite-dropping fireball observations to be detectable. Thus, one likely explanation is that we simply do not have enough observations yet.

Another reason is that even on an orbit with a perihelion distance near the Earth, the impact timescales of the meteoroids produced from these tidal disruptions or other mechanisms (e.g., meteoroid impacts; Turner et al. 2021; Shober et al. 2024a) are just too prolonged. The inverse of typical impact frequencies (10⁸–10¹⁰ yr; Bottke et al. 1994) with the Earth are significantly longer than the decoherence lifetimes of any possible stream. Also, the tidal disruption frequency, estimated to be once every ~2500 years (Schunová et al. 2014), are 10–1000× shorter than the dynamical lifetimes in the inner solar system (Gladman et al. 1997). Thus, if tidal disruption is a consistent mechanism of meteoroid and small asteroid production near the Earth, due to the inherent shortness of stream decoherence, a large majority of the debris being detected telescopically or by fireball networks would necessarily not be expected to be found in distinct streams. This could likely be why Jopek (2020) found many NEOs, generally clusters of q ~ 1 au orbits, an observation we confirm here. Many of these NEOs do not fulfill our statistical significance test; however, the general clusters on these orbits could be representative of older, no longer coherent streams. In addition, the decoherence lifetimes of streams (typically 10⁴-10⁵ years), tend to be significantly shorter than the cosmic-ray exposure (CRE) ages of meteorites (Figs. 15, 16, and 17). Thus, any proposed meteorite linkages in the future must also have a shorter CRE age than the expected decoherence lifetime of a stream on a nearby orbit. This is precisely the reasoning of Shober et al. (2024a), which found that carbonaceous meteorite falls Sutter’s Mill, Flensburg, and Winchcombe have CRE ages as younger or significantly younger than the average CRE age for debris released already on an Earth-crossing orbit.

The only way to have longer CRE ages and still be ejected from an NEO would be to collect much of the CRE age on the surface of the NEO. This would leave a “complex” exposure signature in the irradiation history (Wieler & Graf 2001). One of the primary reasons the Příbram-Neuschwanstein pair was disregarded in the first place was that the CRE ages and meteorite types were significantly different, Příbram being an H5 with a CRE age of ~12 Myr and Neuschwanstein being an EL6 with a CRE age between 43–51 Myr (Stauffer & Urey 1962; Zipfel et al. 2010; Pauls & Gladman 2005). However, as observed with the fall of Almahatta Sitta (asteroid 2008 TC3), this original logic that meteorites originating from the same meteoroid impact necessarily need to share a common composition and irradiation history is understood to no longer be a requirement (Bischoff et al. 2010). Thus, a pair like Příbram-Neuschwanstein could be related and generated during a disruption of a polymict near-Earth rubble pile (Tóth et al. 2011; Granvik & Walsh 2024), assuming a complex exposure is consistent with the cosmo-genic radionuclide and noble gas measurements (Meier et al. 2022). Nevertheless, since the sole reason for the association is orbital similarity, without some additional constraint based on the meteorites themselves, this pair is unfortunately incapable of confirmation. For bodies separated in space for more than a few centuries, orbital similarity discriminants can only detect streams of many objects, not pairs. Even when a stream is well identified, there will always be some sporadic component erroneously included (Shober & Vaubaillon 2024), but without some other constraint, you cannot be sure which objects expressly are the false positives, only the proportion of the sample.

Fig. 16 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 10° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. Decoherence time was defined as when a fictitious meteoroid stream lost 95% or more of the original stream.

Fig. 17 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 20° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. Decoherence time was defined as when a fictitious meteoroid stream lost 95% or more of the original stream.

5 Conclusion

In this study, we revisited the meteoroid and NEO stream deco-herence concept and critically assessed the limits of orbital similarity measures in stream identification. Our findings emphasize the transient nature of near-Earth streams and challenge the traditional reliance on orbital similarity discriminants for individual pair associations. The primary results of this research are summarized as follows:

• Using a Kernel Density Estimation (KDE) method (Shober & Vaubaillon 2024), we evaluated the statistical significance of orbital pairings within various datasets, including 50 meteorite falls with orbits, 616 potential >1 g meteorite-dropping fireballs (of which 350 are likely >50 g), 310 impact detected by US government sensors, and 35 012 NEOs. Our analysis revealed no statistically significant streams in the meteorite fall, USG impact, or fireball datasets. The cumulative D-value distributions for these datasets were consistent with random chance, with no observed clustering beyond what would be expected randomly;

• Conversely, we identified 11 statistically significant clusters within the NEA population and one of 70 fragments associated with the disintegration of comet 73P/Schwassmann-Wachmann. All clusters are found in orbits with perihelion distances near 1au. These clusters, identified using the DBSCAN algorithm with D_H values below 0.03, likely result from tidal disruption events, consistent with the findings of Granvik & Walsh (2024). These clusters suggest that while streams may rapidly decohere, the process generating them is frequent enough to maintain small identifiable NEO families;

• Within this study, we utilized traditional orbit similarity discriminants; however, these orbital discriminates are simplistic and are not ideal for clustering orbital data in all cases. The orbital evolution of streams in the solar system depends on the locations of the streams and the primary forces driving this evolution (which could be size-dependent). The weightings of the orbital parameters within orbital similarity discriminants do not reflect this. Thus, new methods that take into account this orbital evolution variation should be developed to more accurately identify and characterize groupings of comets, asteroids, and meteoroids. There is a statistically significant amount of similarity within the NEO population, but the NEO families identified here using current orbital similarity discriminants should be re-evaluated using a method that takes into account the local dynamics of streams;

• Furthermore, none of the meteorite-dropping fireballs examined (616 possible droppers + 50 meteorite falls) fulfilled the DBSCAN requirements to be members of any of the 12 NEO clusters identified. Even when the minimum cluster size was lowered from five to three members, no impact data was still classified as being member of one of the identified NEO clusters;

• Lyapunov lifetime analysis demonstrated that orbits in near-Earth space, particularly those with semi-major axes between 0.5 and 3.5 AU and eccentricities between 0.2 and 0.9, exhibit chaotic behavior with characteristic lifetimes as short as 60 to 200 years (Fig. 14). These short timescales underscore the rapid divergence of initially similar orbits, rendering long-term orbital similarity measures unreliable for individual pair identification;

• The decoherence lifetimes of meteoroid streams, as defined as the time required for a stream to lose 95% of its original members, ranged from 10⁴ to 10⁵ years for Earth-crossing orbits. This aligns with previous studies such as Pauls & Gladman (2005), indicating that streams on these orbits lose coherence on timescales significantly shorter than most meteorites’ cosmic ray exposure ages. Thus, most meteorites would require the presence of some “complex” cosmic-ray exposure signature to be considered as possibly ejected recently by an immediate-precursor body in near-Earth space.

Our results highlight a significant limitation in using orbital similarity measures: they are not reliable for identifying specific parent bodies of meteorites. The Lyapunov characteristic timescales in near-Earth space are only a few hundred years, after which orbits diverge chaotically. This rapid separation in nearby orbits means that any association between a specific fireball, meteorite, and asteroid, based solely on orbital similarity, is likely coincidental rather than indicative of a true physical connection. Orbital similarity can provide a potential avenue to identify objects with a possible “genetic” link; however, more constraints are always needed to confirm this. While these metrics help identify streams, they cannot confirm specific body associations alone. Despite the evidence that meteorites can be and are ejected from immediate precursor bodies already within near-Earth space (Granvik & Walsh 2024; Shober et al. 2024a), these streams become decoherent rapidly. Identifying such streams using similarity alone within the meteorite and fireball populations will require a significantly larger dataset than we currently have. We would likely need tens of thousands of meteorite-dropping fireball observations to detect these small ephemeral streams with enough statistical significance to separate them from random associations with the sporadic complex if the stream-sporadic ratio is similar to NEOs. Only with such extensive data could we begin to discern small streams resulting from tidal disruptions or other mechanisms, such as meteoroid impacts. In any case, the notion of confidently linking a single fireball to a single NEO through orbital similarity remains fundamentally flawed and unattainable due to the inherent limitations imposed by the chaotic nature of orbits in near-Earth space.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 945298 ParisRegionFP. The Global Fireball Observatory and data pipeline is enabled by the support of the Australian Research Council (DP230100301, LE170100106). FRIPON was initiated by funding from ANR (grant no. 13-BS05-0009-03), carried by the Paris Observatory, Muséum National d’Histoire Naturelle, Paris-Saclay University and Institut Pythéas (LAM-CEREGE). VigieCiel was part of the 65 Millions d’Observateurs project, carried by the Muséum National d’Histoire Naturelle and funded by the French Investissements d’Avenir program. FRIPON data are hosted and processed at Institut Pythéas SIP (Service Informatique Pythéas), and a mirror is hosted at IMCCE (Institut de MécaniqueCéleste et de Calcul des Éphémérides / Paris Observatory). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir superv ised by the Agence Nationale pour la Recherche. This research used Astropy, a community-developed core Python package for Astronomy (Robitaille et al. 2013). Simulations in this paper used the REBOUND software package¹⁰ (Rein & Liu 2012).

References

All Figures

	Fig. 1 Orbital distribution of 50 recovered meteorite falls (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.5. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.
In the text

	Fig. 2 Orbital distribution of 35 012 NEOs (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.05. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.
In the text

	Fig. 3 Orbital distribution of 310 US Government sensor impacts (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.
In the text

Fig. 4 Orbital distribution of 616 possible meteorite falls detected by GFO, FRIPON, or EFN sensors (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15 and ensuring any sample crossed the torus traced out by the Earth’s orbit. The meteorite falls were identified using α − β parameters as described in Sansom et al. (2019) using a minimum end mass of 1 g. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

In the text

Fig. 5 Orbital distribution of 350 probable 50 g meteorite falls detected by GFO, FRIPON, or EFN sensors (blue) versus KDE estimated PDF (red) using a Gaussian kernel with a bandwidth of 0.15 and ensuring any sample crossed the torus traced out by the Earth’s orbit. The meteorite falls were identified using α - β parameters as described in Sansom et al. (2019) using a minimum end mass of 50 g. This PDF, approximating the corresponding sporadic population, was used to estimate the degree of random association in the population via Monte Carlo simulations.

In the text

Fig. 6 Number of meteorite-meteorite pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and DH. The blue circles represent the mean number of pairs from 500 iterations of KDE sampling using a bandwidth of 0.1 based on 350 possible/likely 50 g meteorite falls observed by the GFO, FRIPON, and EFN networks. The blue-shaded area is the 3 σ range possible from random associations. The actual number of low-D value pairs is consistent with the predicted number based on random associations.

In the text

Fig. 7 Number of meteorite-NEO pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H , D′, and DH. The blue circles represent the mean number of pairs identified between 50 random samples of “meteorite falls” with the observed 35012 NEOs from 500 iterations ofKDE sampling using a bandwidth of0.15 forthe350possible50gmeteoritefalls (Fig. 5). Theblue-shaded area is within three standard deviations for random samples of 50 (the number of meteorites with orbits). The actual number of low-D value pairs correctly follows the predicted number based on random associations, indicating no evidence supporting the presence of streams in the recovered meteorite fall dataset. Two asteroids were removed as they correspond to two recovered meteorite falls, Motopi Pan and Almahata Sitta, observed telescopically before impact.

In the text

Fig. 8 Number of CNEOS USG Sensors impact – NEO pairs, ND, with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and DH. The blue circles represent the mean number of pairs identified between 310 random samples drawn from the USG impactors PDF (Fig. 3) with the observed 35 012 NEOs from 500 iterations of KDE sampling. The blue-shaded indicates a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations using the KDE distribution shown in Figs. 2, indicating no substantial evidence supporting streams in the dataset.

In the text

Fig. 9 Number of unique pairs, ND, within the 616 possible 1 g final mass fireballs observed by the GFO-FRIPON-EFN with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for DS H, D′, and Dh. The blue circles represent the mean number of pairs identified between two 616 random samples drawn from the possible 1 g fall PDF (Fig. 4) from 500 iterations of KDE sampling. The blue-shaded indicates a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations, indicating no substantial evidence supporting streams in the fireball dataset.

In the text

Fig. 10 Maximum false positive heat map for bins of semi-major axis (au) by inclination (deg). The heat map coloration represents the mean number of pairs plus 3σ that have DH < 0.1. The sporadic orbital distribution was estimated using the 616 possible 1 g fall observations PDF (Fig. 4). The means and standard deviations were calculated through 500 iterations of KDE sampling. Subplot a examines the maximum number of expected random associations within a population of 616 fireballs, only one (red) in the observations has a measured DH < 0.1, but it is not statistically significant. Subplot b re-examines the significance of the Príbram-Neuschwanstein pair (black square); however, the pair is still not statistically significant within a fireball population of 481 orbits.

In the text

Fig. 11 Number of NEO – NEO pairs with D values below a limiting threshold (red line) compared to the predicted number based on random associations (blue circles with error bars) for Dsh, D′, and Dh. The blue circles represent the mean number of pairs from 500 iterations of KDE sampling using a bandwidth of 0.15, with a shaded blue region indicating a zone encompassing three standard deviations. The number of low-D value pairs corresponds very well to the predicted number based on random associations for large values; however, at very low D values, there appear to be tens of objects with statistically significant similarities. The observational biases of the dataset do not predict the elbow in the cumulative D-value plots.

In the text

Fig. 12 NEO orbital distribution with ephemeris taken from NASA Horizons with clusters identified, objects not belonging to a cluster can be seen as small gray points. The background heatmap, however, indicates the observational bias of the NEO dataset. The red regions denote orbits for which more objects exist in the NEO database relative to the debiased dataset of Granvik et al. (2018), and blue is where there exists a deficit of NEO detections. Clusters identified by Jopek (2020) can be seen as large dark gray points. The clusters identified within this study are colored, each with at least 5 members. They were identified using a DBSCAN algorithm where core points were defined as having at least two associations with a minimum ϵ corresponding to DH value of 0.03. This D value was chosen as it corresponds well to where the cumulative D-value distribution displays a “kink” corresponding to an excess of similarity. The points that were also identified by Jopek (2020) as being in a cluster have a black border. The new large yellow cluster identified is not asteroidal but corresponds to all the fragments of 73P/Schwassmann-Wachmann, which was a Jupiter-family comet observed to have undergone significant fragmentation between 1995 and 2006 (Reach et al. 2009).

In the text

	Fig. 13 Cumulative D-value distributions for every possible pair within 310 USG Sensors orbits 3290 FRIPON fireballs, 824 EFN fireballs, 1913 FRIPON sporadic subset fireballs, 457 EFN sporadic subset fireballs, 50 recoveredmeteoriteorbits, 616potential1 gmeteoritedroppers observed by GFO-FRIPON-EFN, and 35 012 NEOs taken from the NASA Horizons ephemeris service.
In the text

	Fig. 14 Lyapunov characteristic lifetime as a function of the semi-major axis and eccentricity with an inclination of 0° (a), 10° (b), 20° (c). The black lines indicate where the aphelion and perihelion distances equal the nominal semi-major axis value of Mercury, Venus, Earth, Mars, or Jupiter.
In the text

Fig. 15 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 0° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. The ϵ value for the DBSCAN algorithm was chosen to where the “kink” is in the cumulative D-value distributions for NEOs (Fig. 13). Also, core points have at least 2 connections. The decoherence lifetime was thus defined as when the fictitious meteoroid stream lost 95% or more of the original stream mass.

In the text

	Fig. 16 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 10° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. Decoherence time was defined as when a fictitious meteoroid stream lost 95% or more of the original stream.
In the text

	Fig. 17 Decoherence lifetime as a function of the semi-major axis and eccentricity with an inclination of 20° using DS H (a), D′ (b), or DH (c) values in conjunction with a DBSCAN algorithm to determine clustering. Decoherence time was defined as when a fictitious meteoroid stream lost 95% or more of the original stream.
In the text