© The Authors 2026

1 Introduction

The near-Earth asteroid (NEA) (469219) Kamo‘oalewa is the target of the Tianwen-2 sample-return mission (Zhang et al. 2021, 2025) of the China National Space Administration (CNSA), which was successfully launched on 29 May, 2025. Kamo‘oalewa has peculiar characteristics in both its orbit and known physical properties. This object is in a 1:1 mean-motion resonance with the Earth, and it is currently in a quasi-satellite configuration (de la Fuente Marcos & de la Fuente Marcos 2016a). On a timescale of 10 kyr it repeatedly switches between quasi-satellite and horseshoe configuration (de la Fuente Marcos & de la Fuente Marcos 2016a; Fenucci et al. 2022), and its dynamical stability properties makes it a companion of Earth for at least the next 0.5 Myr (Fenucci & Novaković 2021). The available photometric measurements indicate an absolute magnitude of H = 24, meaning that its diameter can be roughly estimated between 30 and 100 m in diameter, and its rotational period is only about 28 minutes (Sharkey et al. 2021). Because of the long observational arc it was possible to determine the Yarkovsky effect (Liu et al. 2022; Fenucci et al. 2025) with a good signal-to-noise ratio. This permitted to also evaluate the thermal inertia with the ASTERIA method (Novaković et al. 2024), indicating a value of ${150}_{-45}^{+90}$ J m⁻² K⁻¹ s^−1/2 (Fenucci et al. 2025). Initial spectral data in the visible and near-infrared wavelengths suggested that Kamo‘oalewa is either an S-type or L-type asteroid (Reddy et al. 2017). Additional measurements in the zJHK colors revealed a high red slope (Sharkey et al. 2021) typical of materials that suffered prolonged exposure to space-weathering effects (Pieters et al. 1993).

The peculiar red slope also resembles that of lunar regolith, which, together with its small size, yields the hypothesis that Kamo‘oalewa originated from the Moon as ejecta from an ancient asteroid impact. Numerical simulations showed that there are rare pathways bringing lunar ejecta into the Earth’s co-orbital space (Castro-Cisneros et al. 2023), which can later transition to a quasi-satellite configuration. By studying the age of the craters on the Moon and their sizes, the Giordano Bruno crater was proposed as the only possible crater from which Kamo‘oalewa could originate (Jiao et al. 2024). Statistics presented in Jiao et al. showed that the impact that created the Giordano Bruno crater could produce 0.3-1 Earth co-orbitals on average. These fragments later need to transition into a quasisatellite state in order to match the dynamical properties of Kamo‘oalewa. Jiao et al. (2024) also compared the available spectral data of Kamo‘oalewa to the Apollo 14 and the Luna 24 mission samples, as well as a representative lunar meteorite. In particular, the Luna 24 sample has been suggested to contain ejecta material from the Giordano Bruno crater (Basilevsky & Head 2012). In the comparison made using the mixture of common factor analyses, the spectrum of Kamo‘oalewa appears to be similar to those of the Apollo 14 samples and the lunar meteorite. On the other hand, the Luna 24 sample is the farthest one from Kamo‘oalewa, which seems in contradiction with the fact that this sample is supposed to contain Giordano Bruno ejecta.

On the other hand, the spectrum of the Giordano Bruno crater was extensively studied by Ogawa et al. (2011) by using data from the Japanese Selenological and Engineering Explorer (SELENE) lunar orbiter. The authors presented the reflectance spectrum between 0.5 and 2 μm wavelengths of several locations across the crater, including the field of ejecta right outside the crater, the rim, the wall, and the floor of the crater itself. All the examined locations show two prominent absorption bands at 1 μm and at 1.3 μm, with significant variations between these two wavelengths depending on the location. Moreover, the steep red slope is not present in any of these data analysed, raising further questions concerning the origin of Kamo‘oalewa.

Laboratory analyses showed that laser-irradiated LL-chondrite powder (size smaller than 45 μm), which simulates the effect of space weathering, can reach a steep spectral slope (Zhang et al. 2024), suggesting that Kamo‘oalewa could be an asteroid that originated in the main belt and suffered from long exposure to space weathering. The same study suggested that LL chondrites, which are thought to originate from the Flora family (Vernazza et al. 2007), are meteorite analogs for Kamo‘oalewa. However, the contribution of the main belt to small NEAs reaching a quasi-satellite state with the Earth has not yet been estimated.

These considerations underscore the need for cautious interpretation of spectral data and that other possibilities should be considered - together with the statistics from the possible dynamical evolution. The upcoming observations of Kamo‘oalewa by the James Webb Space Telescope (JWST) with the Near Infrared Spectrograph (NIRSpec) planned in the first quarter of 2026 (Sharkey et al. 2025) will provide additional high-resolution spectral data, together with better constraints on its size and albedo. The JWST recently demonstrated exceptional capabilities in asteroid studies (Burdanov et al. 2025; Rivkin et al. 2025a,b), and this new dataset is expected to better establish the comparisons with lunar samples ahead of the arrival of the Tianwen-2 spacecraft.

In this work, we estimated the number of Kamo‘oalewa-like objects - i.e., Earth quasi-satellites with a size similar to that of Kamo‘oalewa - originating from the main belt. The relative contribution of ejecta from the Giordano Bruno crater was also evaluated for comparison. We also performed survey simulations to give quantitative estimates of the expected completeness of the population of Kamo‘oalewa-like objects.

2 Methods

2.1 Average fraction of Earth quasi-satellites

We define the Earth’s co-orbital space as orbits with osculating semi-major axis a, eccentricity e, and inclination i, such that $$0.994\text{\hspace{0.17em}au}<a<1.006\text{\hspace{0.17em}au},0<e<0.2,0\text{\hspace{0.17em}deg}<i<20\text{\hspace{0.17em}deg}.$$(1)

On the other hand, we do not put constraints on the longitude of the node Ω, the argument of the pericenter ω, or the mean anomaly M. The resonant dynamics is then identified by the time evolution of the angle σ = λ - λ_E = (Ω + ω + M) - (Ω_E + ω_E + M_E), where the subscript E denotes the corresponding orbital element of the Earth.

We note that, while orbits with larger eccentricities and inclinations could still be Earth co-orbitals, the constraints of Eq. (1) are consistent with the secular evolution of Kamo‘oalewa shown in de la Fuente Marcos & de la Fuente Marcos (2016a) and Fenucci et al. (2022). In fact, over 10 kyr of dynamical evolution, which is a timespan after which the orbit of Kamo‘oalewa becomes stochastic, the semi-major axis of Kamo‘oalewa oscillates between 0.995 au and 1.005 au, while the eccentricity does not exceed 0.14 and the inclination does not grow larger than 10 deg.

To estimate the fraction of Earth co-orbitals that are in a quasi-satellite state, we proceeded with a Monte Carlo approach. We sampled the Earth’s co-orbital space with a large number, N, of initial conditions and numerically propagated them for a time T. All the initial orbital elements are assumed to be uniformly distributed within their domain of definition, since the whole phase space has to be explored. For each initial condition, we determined the time T_QS spent in a quasi-satellite state, and we then computed the average fraction spent in a quasi-satellite state as $$F=\frac{1}{N}\sum _{k=1}^N\frac{T_{QS}^{(k)}}{T}.$$(2)

The number F provides the average fraction of quasi-satellites within the Earth’s co-orbital space at any given time. An averaged approach is needed because Earth co-orbitals may switch their state multiple times, as is also seen in the dynamics of Kamo‘oalewa.

2.2 Numerical simulations and quasi-satellite identification

We fixed the initial epoch at 60 800 MJD, and initial conditions for all the planets from Mercury to Neptune were obtained using the JPL Ephemeris DE441 (Park et al. 2021). Note that the resonant averaged motion is determined only by the initial a and σ (Fenucci et al. 2022; Pan & Gallardo 2025), and, since there are no constraints on ω, Ω, and M, all the possible values of σ were explored via the Monte Carlo approach. Thus, fixing a common starting epoch for the numerical simulations is fully justified by the properties of the resonant averaged dynamics, even if we did not use an averaged dynamical model for the simulations. Numerical integrations were performed with a hybrid symplectic integrator (Chambers 1999), which is able to solve close encounters with planets while maintaining the symplectic characteristics of the integration scheme. This integration method is included in the mercury package (Chambers & Migliorini 1997). While the Yarkovsky effect could be added to the dynamical model of mercury (Fenucci & Novaković 2022), it has been shown that it is not relevant in the short-term dynamical evolution of Kamo‘oalewa (Fenucci & Novaković 2021); hence, we did not include it in the model. We used a time step of 0.5 days for the numerical integrations, and the dynamical state was recorded with a step of one year, which is enough to determine the resonance type of the orbit.

The time spent in a quasi-satellite state, T_QS, was determined using a running-windows approach. First, we computed the time evolution of the resonant argument, σ. To identify a quasi-satellite state, we checked the following conditions on σ and the semi-major axis, a : 1) −50 deg < σ < 50 deg; 2) a must attain values larger and smaller than 1 au at different epochs; 3) σ must have at least two local minima and two local maxima. Using a fixed time window with a length of T_W, all the intervals [t_k, t_k + T_W] were checked for quasisatellite states, where t_k are the output times from the numerical integrations. To catch different libration periods, we scan the whole integration time span with running windows of T_W = 10,20,50,100,150,200,250,300,350 years, and T_W = 400,500,..., 1000 years. The time, T_QS, is then obtained as the sum of the quasi-satellite instances found. Figure 1 shows the time evolution of the resonant angle, σ, for an example orbit, with the green shaded areas indicating the quasi-satellite states identified by the running-windows procedure.

Fig. 1 Time evolution of the resonant angle σ (blue dots) of an example orbit that switches between horseshoe and quasi-satellite states. The green shaded areas are the quasi-satellite states identified with the running-window procedure.

2.3 Expected number of Kamo‘oalewa-like objects

The expected number, N_coorb, of Earth co-orbitals in the same size range as Kamo‘oalewa can be obtained by NEA population models, such as those by Granvik et al. (2018) or Nesvorný et al. (2024a,b). These models provide a size-frequency distribution (SFD) of NEAs, as well as the expected distribution of their orbital elements (a, e, i), thus allowing filtering with the constraints of Eq. (1). We can also assume Poisson statistics on N_coorb; therefore, its standard deviation is $\sigma (N_{\text{coorb}})=\sqrt{N_{\text{coorb}}}$. Note that these NEA population models are based on the migration of asteroids from the main belt to the NEA region by the combined action of resonances and the Yarkovsky effect; thus, all the predicted NEAs are assumed to originate in the main asteroid belt (Granvik et al. 2017). The average number of Kamo‘oalewa-like objects¹ in the steady-state NEA population is thus obtained as $$N_{\text{K/NEA}}=F\times N_{\text{coorb}}.$$(3)

Similarly, we can evaluate the relative contribution of Kamo‘oalewa-like objects originated as ejecta from the Giordano Bruno crater. Jiao et al. (2024) showed that the impact that created the Giordano Bruno crater would produce about 300 fragments of the same size as Kamo‘oalewa, and that 0.0001-1% of them have a chance of becoming Earth co-orbitals at a 3σ confidence level. Thus, the most optimistic estimate for the number of Earth co-orbitals from Giordano Bruno is N_GB = 3. Therefore, the average number of Kamo‘oalewa-like objects from the Giordano Bruno crater is N_K/GB = F × N_GB.

2.4 Survey discovery efficiency of Kamo‘oalewa-like objects

We performed survey simulations to understand the efficiency in discovering Kamo‘oalewa-like objects. We focused on the two most productive currently operational surveys, i.e., PanSTARRS (MPC code F51, Denneau et al. 2013) and the Catalina Sky Survey (CSS, MPC code G96, Fuls et al. 2023). We also simulated the efficiency for the upcoming Vera Rubin Telescope (Ivezić et al. 2019). Given a set of N asteroids that are currently in a quasi-satellite orbit with the Earth, with orbital elements fulfilling Eq. (1) and absolute magnitude 22 < H < 25, we simulated the observational ephemerides from a given survey using the Aegis software (Fenucci et al. 2024) by the European Space Agency (ESA). Ephemerides for G96 were simulated from 1 January, 2005 to 1 October, 2025, while for F51 they were simulated from 1 January, 2011 to 1 October, 2025, thus covering the whole operational period of these two surveys. For Vera Rubin, we computed ephemerides from 1 January, 2026 to 1 January, 2036, which roughly corresponds to its expected operational lifetime.

To simulate the survey strategy to scan the observable sky and to take into account full-moon nights and realistic weathering effects, we used the available telescope pointing data. Historic pointing data for G96 and F51 are available at the Minor Planet Center² (MPC), although those of F51 are available only until 2022. These are real pointings, and they already take into account real observing conditions experienced in the last decades of operations. Simulated pointing data for the next ten years are also available for the Vera Rubin Observatory, which are freely available³. For the simulations done in this work, we used the latest database of pointings called baseline_v5.0.0_10yrs. We also assumed a limiting magnitude of 21.5 for G96, 22.0 for F51, and 24.0 for Vera Rubin, and that every object moving faster than 10 deg day⁻¹ are not detected by automated pipelines (Jedicke et al. 2025). We also simulated the actual visual magnitude of an object observed by a telescope as V + ∆V, where V is the visual magnitude computed from the ephemerides with the HG model using G = 0.15 (Bowell et al. 1989), and ∆V is a factor due to trailing loss. The ∆V term depends on the plane-of-sky velocity of the object and from the characteristics of the telescope, and we computed it as described in Nesvorný et al. (2024a).

Given the ephemerides of an object at a certain epoch, we considered it detectable by the survey if: (1) the simulated visual magnitude V + ∆V was smaller than the limiting magnitude; (2) the corresponding right-ascension (RA) and declination (Dec) coordinates fell within one of the pointings of the survey of that same epoch; (3) the plane-of-sky velocity was smaller than 10 deg day⁻¹. Running the simulations of a sample of Earth quasi-satellites over the operational lifetime of the survey and computing the cumulative distribution of the detectable objects, gives the efficiency in discovering Kamo‘oalewa-like objects.

Fig. 2 Running percentage of the average fraction of quasi-satellites, computed on an increasingly large sample of orbits.

3 Results

3.1 Expected number of Kamo‘oalewa-like objects

We sampled the Earth co-orbital space defined in Eq. (1) with N = 100 000 randomly generated orbits, and their dynamics was propagated for 10 kyr. To verify that the chosen sample is large enough, we computed the average fraction of quasi-satellites on an increasingly larger sample. Figure 2 shows the results of these computations. When the sample is smaller than about 20 000 orbits, the average fraction of quasi-satellites undergoes some fluctuations, while it stabilises for larger samples and approaches a limit when N > 80 000, reaching F = 0.01396 for the complete sample. This proves that our estimate is not dependent on the size of a chosen sample, provided it is large enough. The value obtained means that, on average, only 1.396% of the objects in the Earth’s co-orbital space (Eq. (1)) are quasi-satellite at any given time. To verify that F does not also depend on the total integration time, we extended all the integrations up to 50 kyr. The final value of F is 0.01224, a change of about 10% and only slightly smaller than that obtained by the 10 kyr integrations, thus indicating that the result does not strongly depend on the total integration time span.

As outlined in Sect. 2.3, the number of Kamo‘oalewa-sized Earth co-orbitals originated as ejecta from the Giordano Bruno crater has been estimated to be N_GB = 3 in the best-case scenario; thus, the corresponding average number of Kamo‘oalewa-like objects is N_K/GB = 0.042. As described above, we can evaluate the expected number of Kamo‘oalewa-like objects from NEA population models. Since we are interested in objects that are also similar to Kamo‘oalewa in size⁴, we filtered absolute magnitude values, H, between 22 and 25. With these constraints, the model by Granvik et al. (2018) estimates N_coorb = 88 ± 9.3, while the NEOMOD2 model by Nesvorný et al. (2024a) gives N_coorb = 73 ± 8.5. The diameter of Kamo‘oalewa is not precisely known, but a priori models based on photometry and population models indicate a size between 30 m and 100 m (Fenucci et al. 2025). The NEOMOD3 model by Nesvorný et al. (2024b) also provides the size distribution of the NEA population, and the constraint 30 m < D < 100 m can also be used to extrapolate the number of NEAs in the Earth’s co-orbital region. Figure 3 shows the cumulative SFD of Earth co-orbitals between 25 m and 100 m, extracted from a sample NEA population generated with NEOMOD3. Using the size constraint of Kamo‘oalewa, we obtain N_coorb = 88 ± 9.3, which is consistent with the model by Granvik et al. (2018).

Using the estimate outlined in Sect. 2, we obtain an average number of Kamo‘oalewa-like objects of N_K/NEA = 1.23 ± 0.13, with 22 < H < 25 and with 30 m < D < 100 m, using the Granvik et al. (2018) and the NEOMOD3 models, respectively. We obtain instead N_K/NEA = 1.02 ± 0.12 Kamo‘oalewa-like objects with 22 < H < 25 by assuming the NEOMOD2 model. The right axis of Fig. 3 also shows how the number of Kamo‘oalewa-like objects varies for different values of the diameter using the NEOMOD3 model as input distribution.

Fig. 3 Cumulative size-frequency distribution of NEAs between 30 m and 100 m in the Earth’s co-orbital region, extracted from the NEO-MOD3 NEA population model. Error bars were computed assuming Poisson statistics. The right y-axis scales to the estimated number of Kamo‘oalewa-like objects. Dashed horizontal gray lines highlight the cases with diameters of 30 and 40 m.

3.2 Survey discovery efficiency of Kamo‘oalewa-like objects

To obtain a sample of objects that are currently quasi-satellites of the Earth, we randomly generated orbits with initial Keplerian elements fulfilling Eq. (1) at epoch 60 800 MJD, and with an initial resonant angle, σ, between −50 and 50 deg. The orbits were propagated for 50 years both in the past and in the future, and only those that were quasi-satellites of the Earth for this whole time span were kept. With this procedure, we generated a set of N = 5000 Earth quasi-satellites, which were assigned a random absolute magnitude, H, between 22 and 25. We then simulated the survey efficiency with the settings described in Sect. 2.4.

The left panel of Fig. 4 shows the cumulative fraction of discovered Kamo‘oalewa-like objects by CSS and Pan-STARRS, as a function of the year. CSS shows an almost linear trend in the number of discoveries, ultimately arriving to about 62% after almost 20 years of operations. The biggest limitation of CSS in discovering Kamo‘oalewa-like objects is the limiting visual magnitude, which is often not sensitive enough to spot faint Earth quasi-satellites. On the other hand, the number of discoveries by Pan-STARRS is very steep in the first four years of operations, and the distribution becomes almost shallow after eight years of operations. The final fraction of discovered Kamo‘oalewa-like objects is about 70%. Results for the Vera Rubin observatory are reported in the right panel of Fig. 4. The behavior is similar to that of Pan-STARRS: a large fraction of Kamo‘oalewa-like objects are discovered within the first four years of operations, with few new discoveries in the years after. Our results indicate that Vera Rubin is expected to reach a completeness limit of 90% by 2030, and 92% by 2036. Thus, any undiscovered Kamo‘oalewa-like object would almost certainly be discovered by the Vera Rubin Telescope.

The left panel of Fig. 5 shows the cumulative fraction of discovered quasi-satellites at the end of the simulation time span as a function of the absolute magnitude, H. In the 22 < H < 25 interval, the fraction of discovered objects shows a roughly linear decreasing trend, showing that smaller Earth quasi-satellites are harder to discover. The completeness of CSS at magnitude 22 reaches about 90%, while it drops to 62% at a magnitude of 25. Pan-STARRS is more efficient at smaller diameters, reaching about 70% completeness at a magnitude of 25. Note that, at the same magnitude as Kamo‘oalewa, Pan-STARRS has a discovery efficiency of about 80%, meaning that it would be able to find an extremely large fraction of Kamo‘oalewa-like objects. On the other hand, Vera Rubin has an efficiency of 90% or better through the whole absolute magnitude interval, showing its extreme efficiency in finding Earth quasi-satellites.

Fig. 4 Cumulative fraction of detected Kamo‘oalewa-like quasi-satellites as function of the year. Left: results obtained for the CSS survey (blue curve) and for the PanSTARRS survey (orange curve). Right: expected efficiency of the Vera Rubin Telescope during its ten-year operation period.

Fig. 5 Detected cumulative fraction of Kamo‘oalewa-like objects obtained at the end of the simulations, as a function of the absolute magnitude, H.

4 Discussion

4.1 Implications on the origin of Kamo‘oalewa

The estimates presented above show that, on average, we expect 1.23 ± 0.13 Kamo‘oalewa-like objects at any time. At a magnitude of 24, the expected discovery efficiency of Earth quasisatellites for Pan-STARRS is of 80%, meaning that we expect 0.80 × 1.23 ~ 1 of such known objects. To understand how this estimate fits with the known data, we started from the known Earth co-orbitals and extracted the asteroids that are currently quasi-satellites. To this end, we queried the NEA database of the ESA NEO Coordination Centre⁵ (NEOCC) to obtain all the known NEAs with absolute magnitudes, H, between 22 and 25 and fulfilling the constraints of Eq. (1). The database was queried on 21 November, 2025. A total of eight Earth co-orbitals were obtained (see Table 1), including Kamo‘oalewa. Their orbits were determined with the Aegis system (Fenucci et al. 2024), and, to check their resonant behaviour, we propagated their orbits with the high-fidelity dynamical model for Solar System small bodies implemented in Aegis. Figure 6 shows the evolution of the critical argument, σ, between year 1880 and 2880. Four of these objects are not in resonance with the Earth, and three of them are currently on a horseshoe orbit, leaving Kamo‘oalewa as the only quasi-satellite currently known in its size range. Note that 2015 SO₂ experiences a quasi-satellite configuration longer than 100 years, but starting around year 2350 (see also de la Fuente Marcos & de la Fuente Marcos 2016b), and that Kamo‘oalewa will switch to horseshoe state shortly before that epoch. This gives credibility to the models’ predictions that only one Kamo‘oalewa-like object should exist at any given time. The figure also highlights the period from 2005 to 2026 in gray, which roughly corresponds to the time span of activity of CSS and Pan-STARRS. Note that all the objects were discovered near their passage close to σ = 0 deg, which corresponds to a small difference in mean longitude of the asteroid and the Earth. Thus, the estimates presented here suggest that the current population of Earth quasi-satellites with the same size and orbital properties of Kamo‘oalewa could already be complete, and it can be explained without relying on the lunar-origin hypothesis.

Going down to an absolute magnitude of 26, there are two other known Earth quasi-satellites: 2025 PN₇ and 2023 FW₁₃ (de la Fuente Marcos & de la Fuente Marcos 2025). The linear decreasing trend of the survey efficiency shown in Fig. 5 cannot be extrapolated to fainter magnitudes, because at some point there would necessarily be an exponential decay. For this reason, we took all the orbits generated to obtain the results of Sect. 3.2, assigned them an absolute magnitude of 26, and re-ran the survey simulations for Pan-STARRS. We thus note a completeness limit of 14% reached in 2022. On the other hand, NEOMOD3 predicts a 2427 ± 50 Earth co-orbitals with H < 26, and, since the estimate of Sect. 3.1 does not depend on the size, this converts to about 30 expected Earth quasi-satellites. Combining this number with the discovery efficiency of Pan-STARRS, we would expect 30 × 0.14 ~ 4 known Earth quasi-satellites with H < 26. This is roughly in agreement with the known population, which can again be explained without relying on lunar ejecta.

The results obtained for the Vera Rubin Observatory show that this survey will provide data to assess the validity of our estimates. On the other hand, the general population of Earth coorbitals is more difficult to discover from ground-based observatories alone (see also Morais & Morbidelli 2002), because the σ period of horseshoe and circulating orbits is typically on the order of hundreds of years. The upcoming space-based infrared telescopes NEO Surveyor (Mainzer et al. 2023) by NASA and NEOMIR by ESA (Conversi et al. 2023) may be able to discover a larger part of the remaining Earth co-orbital population, since they are designed to observe at lower solar elongations.

In contrast, the optimistic estimates presented in Sect. 3.1 yield an optimistic number of 0.042 Kamo‘oalewa-like objects originating as Giordano Bruno ejecta, a number that is even slightly lower if we take into account the efficiency of the currently operating surveys. Thus, these numbers suggest that a main-belt origin is statistically more likely than the Giordano Bruno origin hypothesis, and this is consistent with the long space-weathering exposure scenario proposed by Zhang et al. (2024), which therefore appears plausible. Initial in situ observations from the Tianwen-2 spacecrafts should be able to give a better assessment of the origin of Kamo‘oalewa, and we anticipate that a definitive answer will be provided by laboratory analyses of the samples returned to Earth. The samples would allow a better and deeper comparison with the Luna 24 samples containing ejected material from the Giordano Bruno crater (Basilevsky & Head 2012), enabling direct testing of this particular hypothesis on the origin.

If Kamo‘oalewa is actually a lunar ejecta as proposed by Sharkey et al. (2021) and Jiao et al. (2024), it would be the first ever detected NEA to originate from the Moon, marking an extraordinarily low-probability discovery that would open new paradigms for NEA formation and population models (see Sect. 4.2). On the other hand, if Kamo‘oalewa originated in the main belt, as was shown here to be more likely, then the next challenges would be to reconstruct its past evolution, identify possible parent bodies in the main belt, and eventually corroborate the LL-chondrite analog through laboratory analyses of the returned samples. This possibility would also provide a better understanding of the formation process and evolution of small asteroids, which are still poorly constrained, giving a more comprehensive understanding of our Solar System.

Fig. 6 Evolution of the resonant argument, σ, of the eight known Earth co-orbitals with 22 < H < 25 between 1880 and 2880. Each panel shows an angle between −180 and 180 deg. The dashed black line corresponds to σ = 0 deg in each panel. The gray shaded area highlights the period from 1 January 2005 to 1 January 2026.

4.2 Lunar ejecta in the context of NEA population models

Current NEA population models assume that the main belt is the principal source of NEAs. Other low-order effects are also included, such as the catastrophic disruption at a low-perihelion distance (Granvik et al. 2016), and the tidal disruption at close encounters with terrestrial planets (Granvik & Walsh 2024; Nesvorný et al. 2024a). Note, however, that the tidal disruption only has a significant effect at H > 25, i.e., at sizes smaller than that of Kamo‘oalewa. On the other hand, lunar ejecta are not accounted for as a source of NEAs. It is worth noting that the calibration of the model by Granvik et al. (2018) was performed with NEA data obtained between 2005 and 2012, all from the CSS survey. In this time span, the known NEA population comprised fewer than 10 000 NEAs; thus, it was much more limited than that of today. The only asteroids proposed to be lunar ejecta so far, i.e., Kamo‘oalewa, 2020 CD₃ (Bolin et al. 2020) and 2024 PT₅ (Kareta et al. 2025), had not yet been discovered at the time the model was developed. Moreover, Wu et al. (2026) recently estimated the expected annual number of detections of NEAs of lunar origin. The authors obtained 0.75 annual discoveries for Pan-STARRS and 0.31 for the Asteroid Terrestrial-impact Last Alert System (ATLAS). While CSS was not taken into account in their study, it is reasonable to assume a number between 0.31 and 0.75 for CSS. This means that the dataset used to calibrate the model by Granvik et al. (2018) would possibly contain ~5 NEAs of lunar origin, which is probably too few to have introduced any significant bias in the calibration of the NEA population model. The agreement with the NEOMOD3 model calibrated with a larger and more recent sample of data is another hint that lunar ejecta should not constitute a significant fraction of the NEA population, at least at decametres and on larger scales.

Recent estimates by Wu et al. (2026) showed that there could be about 500 000 NEAs larger than 5 m originating from the Moon as impact ejecta, which is still less than 1% of the total expected NEA population. The recently discovered asteroid 2024 YR₄ (Denneau et al. 2024) currently has a probability of 4.3% of impacting the Moon in 2032 (Farnocchia et al. 2026), highlighting that these stochastic events are a current concern (Wiegert et al. 2025). Consequently, although the contribution of lunar ejecta is almost negligible, modeling lunar impact ejecta and incorporating them in NEA population models could be a useful exercise in the near future.

Developments in this direction were recently presented by Jedicke et al. (2025), which put constraints on the distribution of Earth’s mini moons and temporarily captured objects of lunar origin. Although their estimates were affected by significant uncertainties on the frequency of impacts and their parameters, the authors found that lunar ejecta could give rise to ~1-6 temporarily captured objects with a size of ~1 m per year. They also found that this distribution quickly drops as the diameter grows, and the most optimistic estimate for a 20 meter temporarily captured object with lunar origin is one every 1000 years. While this study was focused on Earth’s mini moons and not on the whole Earth co-orbital region, it already shows the quick depletion of the ejecta population at decametres size. For larger asteroids, i.e., larger than ~20 m in diameter, the contribution of lunar ejecta to the whole NEA population has not been clearly quantified. The ejection of fragments of this size typically requires an impactor larger than about 500 meters (Singer et al. 2020). The impact frequency of such objects from the NEA models by Granvik et al. (2018) and Nesvorný et al. (2024b) is ~6 Myr⁻¹ for the Earth. By scaling with the lunar cross-section, the impact frequency for the Moon is ~0.45 Myr⁻¹, which is comparable to the average lifetime of NEAs.

5 Conclusions

Our study provides a quantitative assessment of the potential origin of the asteroid (469219) Kamo‘oalewa through dynamical simulations and current NEA population models. By evaluating the fraction of Earth quasi-satellites within the Earth’s co-orbital region and applying this to main-belt-derived NEA population models, we estimate that the number of Kamo‘oalewa-like objects originating from the main belt is 1.23 ± 0.13, and it is more than one order of magnitude larger than those potentially originating as Giordano Bruno ejecta. We also analysed the discovery efficiency of Kamo‘oalewa-like objects of CSS, PanSTARRS, and the Vera Rubin Observatory. We found a discovery efficiency of between 95% and 70% for an absolute magnitude between 22 and 25 for Pan-STARRS, which is in agreement with the known population of Earth quasi-satellites. The Vera Rubin Observatory should reach a completeness of 95% down to an absolute magnitude of 25, while CSS shows lower discovery efficiency, but it is still greater than 60%.

These findings suggest that Kamo‘oalewa is more likely to have originated from the main asteroid belt, and that it reached the Earth’s co-orbital region later. This interpretation aligns with the recent spectral analyses of laser-irradiated LL-chondrites showing a steep red slope and an absorption band at around 1 μm, both compatible with the data about the spectrum of Kamo‘oalewa. While the lunar-origin hypothesis cannot be entirely ruled out, especially given Kamo‘oalewa’s unique spectral features, our results highlight the importance of a probabilistic and dynamical framework in assessing the origin of NEAs. Observations by the JWST in the first quarter of 2026 are going to provide additional spectral data, and upcoming in situ observations and the eventual sample returned by the Tianwen-2 mission are expected to provide definitive constraints on Kamo‘oalewa’s composition and origin, offering a crucial test of the hypotheses presented here.

Acknowledgements

M.F. dedicates this work to the memory of Professor Àngel Jorba. This research has made use of data and/or services provided by the International Astronomical Union’s Minor Planet Center. P.Z. was supported by the National Natural Science Foundation of China, Grant U24A2008.

References

All Tables

All Figures

	Fig. 1 Time evolution of the resonant angle σ (blue dots) of an example orbit that switches between horseshoe and quasi-satellite states. The green shaded areas are the quasi-satellite states identified with the running-window procedure.
In the text

	Fig. 2 Running percentage of the average fraction of quasi-satellites, computed on an increasingly large sample of orbits.
In the text

	Fig. 3 Cumulative size-frequency distribution of NEAs between 30 m and 100 m in the Earth’s co-orbital region, extracted from the NEO-MOD3 NEA population model. Error bars were computed assuming Poisson statistics. The right y-axis scales to the estimated number of Kamo‘oalewa-like objects. Dashed horizontal gray lines highlight the cases with diameters of 30 and 40 m.
In the text

	Fig. 4 Cumulative fraction of detected Kamo‘oalewa-like quasi-satellites as function of the year. Left: results obtained for the CSS survey (blue curve) and for the PanSTARRS survey (orange curve). Right: expected efficiency of the Vera Rubin Telescope during its ten-year operation period.
In the text

	Fig. 5 Detected cumulative fraction of Kamo‘oalewa-like objects obtained at the end of the simulations, as a function of the absolute magnitude, H.
In the text

	Fig. 6 Evolution of the resonant argument, σ, of the eight known Earth co-orbitals with 22 < H < 25 between 1880 and 2880. Each panel shows an angle between −180 and 180 deg. The dashed black line corresponds to σ = 0 deg in each panel. The gray shaded area highlights the period from 1 January 2005 to 1 January 2026.
In the text