[1]Jose Daniel Castro-Cisneros

[1]Department of Physics, The University of Arizona, Tucson, 85721, AZ, USA

2]Lunar and Planetary Laboratory, The University of Arizona, Tucson, 85721,
AZ, USA

3]Mechanical and Aerospace Engineering, UC San Diego, La Jolla, 92093, CA, USA

# Lunar ejecta origin of near-Earth asteroid Kamo‘oalewa is compatible with
rare orbital pathways

<EMAIL_ADDRESS>Renu Malhotra<EMAIL_ADDRESS>Aaron J.
Rosengren<EMAIL_ADDRESS>*  [  [

###### Abstract

Near-Earth asteroid, Kamo’oalewa (469219), is one of a small number of known
quasi-satellites of Earth; it transitions between quasi-satellite and
horseshoe orbital states on centennial timescales, maintaining this dynamics
over megayears. The similarity of its reflectance spectrum to lunar silicates
and its Earth-like orbit both suggest that it originated from the lunar
surface. Here we carry out numerical simulations of the dynamical evolution of
particles launched from different locations on the lunar surface with a range
of ejection velocities in order to assess the hypothesis that Kamo‘oalewa
originated as a debris-fragment from a meteoroidal impact with the lunar
surface. As these ejecta escape the Earth-Moon environment, they face a
dynamical barrier for entry into Earth’s co-orbital space. However, a small
fraction of launch conditions yields outcomes that are compatible with
Kamo‘oalewa’s orbit. The most favored conditions are launch velocities
slightly above the escape velocity from the trailing lunar hemisphere.

##  1 Introduction

Small bodies in planetary systems can share the orbit of a massive planet in a
long-term stable configuration by librating in the 1:1 mean-motion resonance
(cMsD99, ); such configurations are referred to as co-orbital motion. Examples
of co-orbital arrangements are known for many Solar-System planets, the most
ubiquitous being the large population of Trojan asteroids co-orbiting with
Jupiter. In the context of the idealized circular, restricted three-body
problem (CR3BP), there are three main types of co-orbital states:
Trojan/tadpole (T), horseshoe (HS), and retrograde satellite/quasi-satellite
(QS) (fNetal99, ). The two cases of interest, horseshoe and quasi-satellite,
are shown in Fig. 1 (a), which are distinguished by the center of oscillation
of their longitudes relative to Earth, of $180^{\circ}$ and $0^{\circ}$,
respectively.

Unlike the long-term stable population of Trojan asteroids co-orbiting with
Jupiter, most near-Earth asteroids (NEAs) have chaotic orbits with dynamical
lifetimes much shorter than the age of the Solar System (Gladman:2000, ), and
asteroids stably co-orbiting with the Earth on such timescales are uncommon.
An assessment of Earth’s co-orbital companions shows a total population of at
least twenty-one objects, with two Trojan-type, six in the QS state, and
thirteen undergoing HS motion; all of these objects are in their co-orbital
states only temporarily, typically on less than decadal timescales (mKsC20, ;
cFM21, ; sD23, ). The recently discovered quasi-satellite of the Earth,
(469219) Kamo‘oalewa, is exceptional among the Earth’s co-orbitals due to the
longer-term persistence of its HS–QS transitions (cFM14, ; bD19, ; mReS20, ;
mFbN21, ).

Kamo‘oalewa’s diameter is estimated to be 46–58 m (bs21, ), and its orbital
elements are listed in Table 1, in which we observe that, although its semi-
major axis is very close to Earth’s, its orbital eccentricity and inclination
are not atypical of NEAs. Its ephemeris over a few centuries in the past and
in the future, obtained from the Jet Propulsion Laboratory’s (JPL) Horizons
web service, shows that the transition from HS to its current QS state
occurred nearly a century ago; an event that will reverse in about 300 years
when Kamo‘oalewa will again pass into a HS orbit (Fig. 1 (c)).

Long-term numerical simulations indicate that these transitions will recur
over hundreds of thousands or even millions of years (cFM16a, ; mReS20, ;
mFbN21, ). This can be contrasted with Earth’s first-known recurrent quasi-
satellite, asteroid 2002 AA29, whose future predicted QS state will last only
for a few decades (rB04, ; pW09, ). Kamo‘oalewa’s close proximity to Earth and
its unknown dynamical origin make it a scientifically compelling candidate for
a future space mission (cV19, ; wJ20, ).

Figure 1:  Co-orbital dynamics in the three-body problem and its relation to Kamo‘oalewa’s orbital dynamics. (a) The two classes of co-orbitals considered in this work: horseshoe companions oscillate about the $L_{3}$ Lagrange point, diametrically opposite the planet’s location, and encompass both $L_{4}$ and $L_{5}$ Lagrange points; and quasi-satellites orbit outside the planet’s Hill sphere and enclose both the collinear $L_{1}$ and $L_{2}$ Lagrange points. (b) The trace of asteroid (469219) Kamo‘oalewa’s path in Cartesian coordinates in the co-rotating frame; Earth’s position is shown in blue. (c) Kamo‘oalewa’s semi-major axis $a$ and relative mean longitude $\Delta\lambda=\lambda-\lambda_{\mathrm{Earth}}$ as a function of time, with Horseshoe motion appearing in violet, while quasi-satellite motion is shown in green. The orbital propagation data for 1600–2500 CE are from JPL’s Horizons ephemeris service (retrieved on 8 June 2023). Orbital elements |  | Value | Uncertainty
---|---|---|---
Semi-major axis |  $a$ (AU) | 0.9989754217067754 | $3.5408\times 10^{-9}$
Eccentricity | $e$ | 0.1064616822197207 | $2.4405\times 10^{-7}$
Inclination | $i\,(^{\circ})$ | 7.737141555926749 | $1.6932\times 10^{-5}$
Longitude of ascending node | $\Omega\,(^{\circ})$ | 67.69308146089658 | $1.4631\times 10^{-5}$
Argument of perihelion | $\omega\,(^{\circ})$ | 311.1680143115627 | $2.3126\times 10^{-5}$
Mean anomaly | $M\,(^{\circ})$ | 74.35927547581252 | $2.2534\times 10^{-5}$
Table 1:  Orbital elements for Kamo‘oalewa. The orbit has been determined at
epoch J2452996 (22 December 2003) from the JPL Horizons web-service (retrieved
on 8 June 2023).

Several hypotheses have been proposed for the origin of Kamo‘oalewa (cFM16a, ;
bs21, ). Sharkey and colleagues measured its reflectance spectrum and found it
to have an L-type profile resembling lunar silicates bs21 , inconsistent with
typical NEAs. These authors also concluded that Kamo‘oalewa is unlikely to be
an artificial remnant from an earlier lunar mission. Its modest inclination
could be indicative of a temporarily captured NEA, as is speculated for other
planetary co-orbitals (cFM12, ). Its orbital eccentricity, however, is
atypical of such captured co-orbital states found in numerical simulations,
which generally range between 0.3 (Venus crossing) and 0.6 (Mars crossing)
(mMaM02, ). The other proposed scenarios are that Kamo‘oalewa might have
originated in the Earth-Moon system, either from a hitherto undiscovered
quasi-stable population of Earth’s Trojans or as a lunar ejecta from a
meteoroidal impact (bs21, ).

These latter scenarios for the provenance of Kamo’oalewa are at variance with
prevailing theoretical models of near-Earth objects (NEOMOD2018, ; NEOMOD2023,
) as these models assume only the main asteroid belt and comets as sources of
NEAs. As a check, we employed the NEOMOD simulator (NEOMOD2023, ) and found
that the latest model of NEAs does not account for Kamo’oalewa-like orbits.

The focus of the present work is to examine the hypothesis that Kamo‘oalewa
originated as lunar ejecta. We approach this by numerically simulating test
particles (TPs) launched from the Moon’s surface and following their
subsequent orbital evolution. We use a physically plausible range of launching
speeds and directions and four representative launch locations (see Fig. 2).
The dynamical evolution of lunar ejecta has been previously investigated with
numerical simulations bG95 . While those authors focused on determining
whether lunar ejecta impact the Earth or Moon or escape into heliocentric
orbits, our work focuses on determining whether such particles have dynamical
pathways that lead to co-orbital states. This is a more delicate question
because, as we will see, such outcomes require statistically rare initial
conditions (ICs); to our knowledge, this question has not been previously
investigated.

Figure 2: Initial conditions of lunar ejecta. Launch conditions for test
particles in terms of the parameters $\theta_{1}$, $\theta_{2}$, and
${v}_{\text{L}}$. The unit vector $\hat{r}$ defines the direction pointing
towards Earth and $\hat{t}$ the transversal direction.

##  2 Results

As in previous numerical investigations of lunar ejecta, a variety of
dynamical behaviors were found as particles entered heliocentric orbits. In
order to depict the global results graphically, we projected the orbital
evolution of the particles onto the semi-major axis–eccentricity plane, as
shown in Fig. 3.

An immediately observable feature in this diagram is that most of the launched
particles evolve into orbital parameter regions traditionally demarcated as
the Aten and Apollo regimes of the population of near-Earth objects. A similar
result has been reported previously (wB96, ); it supports the suggestion that
some of the members of the Aten and Apollo dynamical groups originate as lunar
debris (cFM16a, ). The other noteworthy feature in Fig. 3 is the vertical
structure of a low density of points around $a=1$ AU. This is the co-orbital
region where we find the current orbit of Kamo‘oalewa and other HS–QS co-
orbital NEAs. The evident well-defined boundaries of this region show that
there is a dynamical barrier between lunar ejecta and co-orbital states but
the finding of some outcomes in this region indicates that the barrier is
somewhat porous, allowing a small fraction of lunar ejecta to evolve into and
remain in co-orbital states for varying periods of time.

Table 2:  Summary of fates of lunar-ejecta TPs. Frequencies of the orbital fates (including collisions with Earth and Moon) of the TPs from the four representative launching sites. Percentages of the total number of launched particles per site is shown in parenthesis. (Note that a TP may have reached a co-orbital state before colliding).  Launch site | # HS | # HS-QS | # Moon colliders | # Earth colliders
---|---|---|---|---
Near side | 189 (5.11 %) | 24 (0.65 %) | 31 (0.84%) | 301 (8.14 %)
Trailing side | 281 (7.60 %) | 42 (1.13 %) | 41 (1.11 %) | 373 (10.10 %)
Far side | 157 (4.24 %) | 24 (0.65 %) | 36 (0.97%) | 303 (8.19 %)
Leading side | 236 (6.37 %) | 31 (0.84 %) | 9 (0.24%) | 107 (2.89 %)

We identified the co-orbital outcomes by visual inspection of the time
evolution of the particles that spend some time within the semi-major axis
zone of 0.98–1.02 au. Overall, we found that 6.6% of all launched particles
exhibited at least temporary co-orbital motion, most as HS (5.8%) and some as
HS–QS (0.8%). A particle had to perform at least one HS or QS oscillation to
be considered temporarily in a co-orbital state. A quantitative summary of the
frequency for each dynamical outcome from each of the four launch sites is
given in Table 2 along with the total collisions detected from each site. The
trailing side produced the most co-orbiters (both HS and QS), followed by
leading side, and next by near side and far side which produced similar
statistics.

Figure 3:  Summary of orbital outcomes of lunar ejecta particles in orbital
parameter space. (a) Orbital evolution for 5,000 years of 14,800 lunar ejecta
particles, projected on the $(a,e)$ plane. (b) The evolution of Kamo‘oalewa
and four different Kamo‘oalewa-like (KL’s) cases are highlighted on a zoomed-
in portion of the diagram. The cadence of the non-co-orbital trajectories is
downsampled (one point per 250 years) for legibility of the plots.

Amongst all the initial conditions we simulated, some are more favorable for
co-orbital outcomes than others. This is illustrated in Fig. 4; the histogram
in Fig. 4 (b) shows the distribution of the initial launch speed of the cases
that resulted in HS and HS–QS outcomes. Overall, the most favored launch
velocity for HS outcomes is near the minimum of the sampled range (i.e., just
above lunar escape velocity); for HS–QS outcomes (i.e., Kamo’oalewa-like) the
most favored initial launch speed is moderately higher, $\sim(4.0-4.4$) km
s-1, in agreement with the estimates from Section 5.1. In general, the total
number of HS–QS outcomes decreases as the speed increases, discouraging
exploration for larger values.

Figure 4:  Outcomes of lunar ejecta particles related to their launch
conditions. (a),(b) Each point represents a launch condition: one of four
launch locations (near-side, far-side, leading side, and trailing side), and
the launch velocity (decomposed into its radial and tangential components).
See Fig. 2 for an illustration of the launch locations and launch velocity
component directions. The points in red highlight those launch conditions that
result in a HS state while the points in dark blue correspond to detected
HS–QS transitions during the 5,000 years of simulation time. (c) Histogram of
the frequencies of co-orbital outputs with respect to the launching speed.

For additional detail, we examined the outcomes as a function of the radial
and tangential components of the launch velocity (see Fig. 2 for an
illustration of those directions), and made scatter plots of these velocity
components for the four launch locations. In Fig. 4 (a), the initial launch
velocity components from the lunar near-side are plotted as the yellow dots
and those from the lunar far-side are in gray. In Fig. 4 (b) , the initial
launch velocity components from the lunar trailing-side are plotted as the
yellow dots and those from the lunar leading-side are in gray. The HS and
HS–QS outcomes are highlighted as the red and blue points. A clear asymmetry
can be observed in these diagrams: most of the co-orbital outcomes were
launched with a negative tangential velocity (i.e., in the trailing direction
of the Moon’s orbital motion). It is also evident that, out of the four
representative launch sites considered, the trailing side is the most prolific
in producing co-orbiters. Additionally it can be noticed that most of the co-
orbitals produced from the leading side arise from the lower launch speeds,
while for the other sites most of them arise from moderately larger launch
speeds ($>\sim 3$ km s-1) (this will be shown to be due to the higher
frequency of collisions for these conditions, as exposed in Fig. 5 ). We can
also observe that for the larger launch speeds, in the range $4$–$6$ km s-1,
co-orbital outcomes are favored for launch directions in the radial or anti-
radial direction. These patterns in the outcomes are consistent with the
theoretical expectations outlined in §5.1.

It is perhaps noteworthy that we did not find any tadpole-type outcomes, that
is, particles librating around just the L4 or the L5 Lagrange points. Other
possible fates that were examined were collisions. Collisions with the Moon
and the Earth were registered, most of them occurring at the lower launch
speeds and within the first 100 years of their evolution. The statistics of
the collision outcomes is shown in Fig. 5, in a scheme analogous to Fig. 4.
That is, panels (a) and (b) show the scatter plots of the initial conditions
that end in collisions and panel (c) plots a histogram of the frequency of
collisions at different launch speeds. We observe a clear, rapid decay of
collision outcomes for larger launch speeds. The distribution appearing in
Fig. 5 (a) and (b), being concentrated at the lower speeds from near, far, and
trailing side, accounts for the low frequency of co-orbital outcomes under
these conditions, as a particle that may have reached such a state would have
collided before it could enter into a co-orbital state.

Figure 5:  Collisional outcomes of lunar ejecta particles related to their
launch conditions. (a),(b) Each point represents a launch condition: the
radial and tangential components of the launch velocity and one of four
representative locations on the lunar surface (near-side, far-side, trailing-
side and leading-side). Points in red indicate collisions with the Moon,
points in blue indicate collisions with Earth. (c) Histogram of the
frequencies of collisions with respect to the launch speed.

Among the cases of HS-QS co-orbital outcomes observed in the simulations, most
of them (around 66%) displayed only one transition or departed the QS state
rapidly (before 1000 years), performing only one or two transitions. The
orbits of interest are those whose HS–QS transitions recur persistently in a
stable fashion for thousands of years, like Kamo‘oalewa. For the nine ICs that
showcased such Kamo‘oalewa-like dynamics (henceforth referred to as KL’s; see
Fig. 3), the evolution was tracked further, for up to 100,000 years, or until
they departed their co-orbital states. Fig. 6 shows an example KL outcome with
persistent HS–QS transitions; this particle has recurrent residence times of
$400$–$600$ years in the QS state, in-between shorter residence times in the
HS state. For comparison, Kamo‘oalewa’s current time of residence in the QS
state is approximately 400 years.

Figure 6: An example of a Kamo’oalewa-like outcome of a lunar ejecta particle.
(a) Evolution of the semi-major axis and (b) relative mean longitude for the
KL1 trajectory (${v}_{\text{L}}=5.1$ km s-1). HS motion is shown in blue,
while QS motion is shown in yellow. The HS-QS oscillations persist for up to
4500 years (not shown).

As previously noted, Kamo‘oalewa possesses a modest ecliptic inclination of
about $8^{\circ}$. The orbital planes of most simulated ejecta particles were
found to remain close to the ecliptic plane (typical inclinations $\sim$
$1^{\circ}$–3∘). This is because of our adopted simplification of considering
initial launch conditions in a projected plane close to the ecliptic (see Fig.
2 and the associated description in §5.2). However, we did find some cases of
KL outcomes in which the inclinations were driven up to higher values,
reaching inclinations similar to Kamo‘oalewa’s. This can be seen in Fig. 7,
where we plot the time evolution of the inclination of a simulated particle
(KL2), as well as that of Kamo‘oalewa’s orbit. We observe episodic higher
inclination states that persist for a few thousand years in-between short-
lived low inclination states. The higher inclination states are reached by a
sequence of inclination jumps that occur at close approaches between the
particle and Earth, and these jumps build up coherently over time spans of a
few hundred years. In this figure, we also plot Kamo’oalewa’s inclination
evolution over a 10,000 year time span, revealing similar inclination jumps at
close approaches with Earth during its HS state as well as similar features
when transitioning from HS to QS states, but with a boost in inclination
rather than a decrease. These results demonstrate that Kamo’oalewa’s
inclination could have arisen from a smaller initial inclination by means of
kicks at close approaches during its HS state.

Figure 7: Evolution of the ecliptic inclination of Kamo’oalewa and of a lunar
ejecta particle. For Kamo’oalewa, HS states appear in violet and QS states in
green. For the ejecta (KL2, ${v}_{\text{L}}=4.4$ km s-1), non-co-orbital
motion is shown in black, HS states appear in blue and QS states in yellow.
The orbital evolution was obtained with numerical propagation under the same
model, starting at initial epoch J2452996.

##  3 Discussion

It has been suggested that three small NEAs — 2020 PN1, 2020 PP1, and 2020
KZ2, of estimated sizes in the range 10–50 m 2020PN1 ; 2020PP1 ; 2020KZ2  —
may have the same provenance as Kamo‘oalewa due to their close orbital
clustering and the similarities they exhibit in their orbital evolutions on
timescales of a few thousand years cFM21 . We have not investigated the
orbital dynamics of these individual objects, but their resemblance to
Kamo‘oalewa’s orbital elements implies that our results for Kamo‘oalewa could
also be applicable to these objects’ origin.

The lunar ejecta hypothesis for the provenance of Kamo‘oalewa and other small
Earth co-orbitals can be tested for consistency with the lunar impact crater
record and cratering mechanics. The lunar ejecta velocities (in excess of
lunar escape speed, 2.4 km s-1) needed to obtain the co-orbital outcomes
appear to be achievable in meteoroidal impacts on the Moon. Impacts on the
Moon have typical impact speed of 22 km s-1 and as high as 55 km s-1
(Ito:2010, ; Gallant:2009, ). Very small ejected debris particles may achieve
comparable speeds, although the total fraction of such very high–velocity
ejecta (solid or molten) is exceedingly small (Liu22, ; Melosh87, ). Based on
studies of lunar secondary craters, it is estimated that an escaping lunar
ejecta fragment of size in the tens of meters would be expected only from
relatively large impact craters, of diameter exceeding $\sim$ 30 km
(Singer:2020, ; M20, ). During the past $\sim$ 1100 Ma of lunar history (the
Copernican period in the lunar geological timescale), there were 44 impact
craters of diameter exceeding 30 km (Stoffler:2006, ), indicating that such
large impacts occur at average intervals of about 25 Myr. The implication is
that if Kamo‘oalewa is a lunar impact ejecta fragment, then it was launched
from the lunar surface ${\cal O}(10^{7})$ years ago. We leave to a separate
study to investigate whether a lunar crater of appropriate size and age and
geographic location can be consistent with the lunar ejecta hypothesis for the
provenance of Kamo‘oalewa. If supported by such studies, Kamo‘oalewa would, to
the best of our knowledge, be the first near-Earth asteroid to be recognized
as a fragment of the Moon. It would be of great interest for cosmochemical
study as a sample of ancient lunar material. The rarity of Kamo’oalewa-like
orbital outcomes (compared to Aten- and Apollo-like outcomes) in our
simulations of escaping lunar ejecta suggests that many other lunar ejecta
remain to be identified amongst the background population of near-Earth
asteroids. This prediction is testable with near-infrared reflectance spectra
of very large numbers of NEAs that will be obtained by the forthcoming Near-
Earth Object Surveyor project (neosurveyor, ).

Additional exploration of the orbital evolution of lunar ejecta is also
warranted. Our numerical investigations reported here were limited in a number
of ways, so it is useful to list some future directions of investigation. In
the present study, we identified the most-favorable launch velocities of lunar
ejecta for Kamo‘oalewa-like outcomes for initial conditions of the Solar
System taken near the present epoch. Although in §5.2 we invoked the
Copernican principle that the present epoch is not “special,” we do recognize
that our results may have some sensitivity to the initial epoch. The most
important limitation is due to Earth’s orbital eccentricity, which is time
variable and undergoes excursions up to five times its current value on
timescales of ${\cal O}(10^{6})$ years. Consequently, Earth’s orbital velocity
varies by $\sim 2$ km s-1, an amount that is comparable to (or a significant
fraction of) the launch velocities of escaping lunar ejecta. Therefore, this
could influence the frequency of the co-orbital outcomes of escaping lunar
ejecta. Thus, sampling initial epochs when Earth’s eccentricity is different
is needed to understand more comprehensively the statistics of co-orbital
outcomes of lunar ejecta. Sensitivity to epoch could also arise from lunar
phase at launch (because the relative magnitude of solar perturbations on
escaping lunar ejecta particles at full moon versus new moon phase is also a
significant fraction of the lunar orbital velocity) and from the perturbations
of Jupiter and other planets that would be slightly different at different
epochs.

We would also like to understand the dynamical mechanism by which
Kamo‘oalewa’s persistent HS–QS transitions occur. Of the three possible co-
orbital states, the QS state is the rarest found among small bodies in the
Solar System. In the simple model of the planar, circular, restricted three
body problem (PCR3BP), the intrinsic stability of nearly coplanar QS orbits
has long been established (ySsT08, ; aP17, ; gVkA18, ). It has been linked to
the existence of the family f of periodic orbits in the PCR3BP (mH69, ), and
referred to as distant retrograde orbits (DROs) in applications to spacecraft
navigation and mission design (cSdS10, ; mL18, ). In the spatial problem,
vertical instabilities can arise and transitions between co-orbital states are
possible (sM06, ; vS14, ; hMfN18, ; kO18, ; yQdQ22b, ). In the regime of large
eccentricity and inclination, (fNetal99, ) has attributed such transitions to
a secular drift of the asteroid’s perihelion and (cFM16a, ) has suggested this
as a mechanism that applies to Kamo‘oalewa. While the eventual escape from co-
orbital states may be linked to planetary secular perturbations (pW10, ;
aCdA11, ; mC14, ), or to planetary close encounters (aCnG21, ), or to
Yarkovsky-driven migration (mFbN21, ), it is likely that the short-time
transport dynamics of Kamo’oalewa are governed by the invariant-manifold
structure of the Lagrange points (eBmO08, ; kO18, ; yQdQ22a, ). For example,
some authors attribute the entry and escape mechanisms of Kamo‘oalewa’s HS–QS
transitions to such phase-space structures (yQdQ22a, ), but others invoke
chaotic tangles of the Lagrange points to explain the dynamical mechanisms of
capture into sticky QS orbits (cSdS10, ). Nevertheless, it is challenging to
identify the specific phase-space structures responsible for the dynamical
transport phenomena exhibited by Kamo‘oalewa-like objects. More research is
needed to understand the precise role of these manifolds on the dynamics of
co-orbital objects like Kamo‘oalewa, as well as on the wider NEA populations,
and their implications for the asteroid impact hazard on our home planet
(nT20, ; yQdQ22b, ).

The complex and non-linear nature of the calculations performed leads to a
large sensitivity to several conditions. For instance, initial conditions for
objects in the Solar System were gathered from the JPL Horizons service, where
masses and orbital elements are subject to updates and refinements. Further
inconveniences arise from the fact that the orbital fates are classified based
on visual inspection of 5000 year of orbital evolution of more than 10,000
simulated particles, so results are vulnerable to human error. Different
results may also arise if the initial integration time of 5,000 years is
modified.

##  4 Conclusions

In our numerical simulations of the dynamical fates of lunar ejecta, we
explored a representative range of ejecta launch conditions expected from
large meteoroid impact events. The vast majority (more than 93%) of the launch
conditions we considered resulted in ejecta reaching heliocentric orbits
similar to the Aten and Apollo groups of NEAs, with no co-orbital behavior
detected; this is consistent with previous results (wB96, ; bG95, ). However,
in a small minority (6.6%) of cases we detected the existence of pathways
leading to co-orbital states, most commonly horseshoe orbits, but also those
resembling Kamo’oalewa’s; the latter exhibit persistent transitions between
quasi-satellite and horseshoe orbits. These minority outcome events have not
been previously reported. The existence of these outcomes lends credence to
the hypothesis that Kamo‘oalewa could indeed be lunar ejecta. The launch
conditions most favored for such an outcome are those with launch velocities
slightly above the lunar escape velocity and launch locations from the Moon’s
trailing side. We also find that Kamo’oalewa’s inclination may have been
boosted by close approaches with the Earth during its horseshoe state.

##  5 Methodology

###  5.1 Theoretical estimates

We begin with the observation that particles originating in the Earth-Moon
(EM) system that escape and evolve into Earth-like heliocentric orbits,
including co-orbital states such as horseshoe and quasi-satellite types, would
be those that escape with low relative velocity with respect to the EM
barycenter. Here we make some estimates of the dynamical conditions of launch
from the lunar surface that would favor outcomes with Earth-like heliocentric
orbits.

For these estimates, we will take the Earth’s Hill sphere as the approximate
boundary between geocentric and heliocentric space, and the lunar Hill sphere
as the approximate boundary between selenocentric and geocentric space. The
radius of the lunar Hill sphere is approximately 35 lunar radii:

$r_{\text{H},{\rightmoon}}=\Big{[}\frac{m_{\rightmoon}}{3\tilde{m}_{\oplus}}\Big{]}^{\frac{1}{3}}a_{\rightmoon}\simeq
35.2~{}R_{\rightmoon},$ (1)

and Earth’s Hill sphere radius is approximately 1% of Earth’s heliocentric
orbit radius:

$r_{\text{H},\oplus}=\Big{[}\frac{\tilde{m}_{\oplus}}{3(m_{\odot}+\tilde{m}_{\oplus})}\Big{]}^{\frac{1}{3}}a_{\oplus}\simeq
0.01~{}\hbox{au}.$ (2)

Here $m_{\rightmoon},\tilde{m}_{\oplus},m_{\odot}$ are the lunar mass, the
mass of Earth + Moon, and the solar mass, respectively,
$a_{\text{\rightmoon}},a_{\oplus}$ are the lunar orbit radius and Earth’s
heliocentric orbit radius, respectively (both approximated as circular
orbits), and $R_{\rightmoon}$ is the lunar radius.

We denote with $\boldsymbol{v}_{\text{L}}$ the launch velocity of a particle
launched from the Moon’s surface; this is relative to the lunar barycenter.
Particles launched with $v_{\text{L}}$ exceeding the Moon’s escape speed,
$v_{\text{\rightmoon,esc}}=2.4$ km s-1, will reach the lunar Hill sphere
boundary with a residual speed $\boldsymbol{\delta v_{\text{L}}}$ relative to
the lunar barycenter. The magnitude of this residual velocity is estimated
from the equation for conservation of energy in the lunar gravitational field,
and is given by

$(\delta
v_{\text{L}})^{2}=v_{\text{L}}^{2}-\frac{2Gm_{\rightmoon}}{R_{\rightmoon}}\Big{(}1-\frac{R_{\rightmoon}}{r_{\text{\text{H}},{\rightmoon}}}\Big{)}=v_{\text{L}}^{2}-0.97v_{{\rightmoon},\text{esc}}^{2}$
(3)

For later reference, we observe that for a particle launched with
$v_{\text{L}}\approx v_{\text{\rightmoon,esc}}$, the residual velocity at the
lunar Hill sphere radius is $\delta v_{\text{L}}\approx 0.4$ km s-1.

The velocity of an escaping lunar ejecta particle relative to the Earth-Moon
barycenter, $\boldsymbol{v}_{\text{EM}}$, is found by adding (vectorially) the
residual velocity $\boldsymbol{\delta v}_{\text{L}}$ to the lunar orbital
velocity, $\boldsymbol{v}_{\rightmoon,\text{orb}}$,

$\boldsymbol{v}_{\text{EM}}=\boldsymbol{\delta
v}_{\text{L}}+\boldsymbol{v}_{\rightmoon,\text{orb}}$ (4)

The Moon’s orbital velocity about the Earth-Moon barycenter is
$v_{\text{\rightmoon,\text{orb}}}=1.0$ km s-1. To escape from the Earth-Moon
system, the magnitude of $v_{\text{EM}}$ must exceed
$\sqrt{2}v_{\text{\rightmoon,\text{orb}}}\simeq 1.4$ km s-1. From Eq. 4, we
see that this requires that $\delta v_{\text{L}}$ must exceed 0.4 km s-1. This
minimum value is coincidentally the same as the residual velocity at the lunar
Hill sphere of lunar ejecta launched with just-the lunar escape speed, that
is, $v_{\text{L}}\approx v_{\text{\rightmoon,esc}}$.

The magnitude of $\boldsymbol{v}_{\text{EM}}$ depends upon the location and
speed of launch. We illustrate with two limiting cases. First consider lunar
ejecta launched in the vertical direction from the apex of the lunar leading
hemisphere. Such ejecta will reach the lunar Hill sphere boundary with
residual velocity in a direction nearly parallel to the Moon’s orbital
velocity. Consequently, they will get boosted by $\sim 1$ km s-1 (the lunar
orbital velocity) to $v_{\text{EM}}>1.4$ km s-1, assuring escape from Earth’s
Hill sphere. The second limiting case is that of lunar ejecta launched
vertically from the apex of the trailing lunar hemisphere. Such ejecta will
reach the lunar Hill sphere boundary with residual velocity approximately
anti-parallel to the Moon’s orbital velocity, so their $v_{\text{EM}}$ will be
lower than $\delta v_{\text{L}}$ by $\sim 1$ km s-1. In this case a residual
velocity of magnitude $\delta v_{\text{L}}>2.4$ km s-1 is needed in order to
achieve $v_{\text{EM}}>1.4$ km s-1. From Eq. 3, this requirement implies a
launch velocity $v_{\text{L}}>3.4$ km s-1. Ejecta from other locations and
different launch directions will require a minimum launch speed in the range
2.4 to 3.4 km s-1 in order to leave the Earth-Moon Hill sphere.

In the geocentric phase, the initial location of an escaping lunar ejecta
particle is approximately at one lunar orbit distance,
$a_{\text{\rightmoon}}$, and its velocity is $v_{\text{EM}}$ relative to the
Earth-Moon barycenter. The ejecta will reach the Earth-Moon Hill sphere
boundary with a residual velocity, $\delta v_{\text{EM}}$, relative to the
Earth-Moon barycenter given by the energy conservation equation in geocentric
space,

$(\delta
v_{\text{EM}})^{2}=v_{\text{EM}}^{2}-2\frac{G\tilde{m}_{\oplus}}{a_{\text{\rightmoon}}}\Big{(}1-\frac{a_{\text{\rightmoon}}}{r_{\text{H},\oplus}}\Big{)}\simeq
v_{\text{EM}}^{2}-(1.2~{}\mathrm{km~{}s}^{-1})^{2}.$ (5)

Taking $v_{\text{EM}}=1.4$ km s-1 (the minimum required to achieve escape from
geocentric space), escaping lunar ejecta will enter heliocentric space with a
residual velocity (relative to the Earth-Moon barycenter) of $\delta
v_{\text{EM}}=0.7$ km s-1. This is a small fraction, $\sim 0.023$, of Earth’s
heliocentric orbital velocity. Particles having $\delta v_{\text{EM}}$ close
to this minimum value will enter heliocentric space with the most Earth-like
orbits, and would be good candidates for entering co-orbital states such as
the horse-shoe or quasi-satellite orbits.

For high speed lunar ejecta, those launched in a direction opposite to the
lunar orbital velocity achieve lower residual velocities relative to the
Earth-Moon barycenter. This means that launch locations from the lunar
trailing hemisphere (that is, the hemisphere opposite to the Moon’s orbital
motion around Earth) would be more favorable for Earth-co-orbital outcomes of
escaping lunar ejecta. The circumstance that the Moon is in synchronous
rotation with its orbital motion (and likely has been so for most of its
history (Wieczorek:2009, ; Aharonson:2012, )) means that the favorable launch
location can be geographically constrained in this way for even ancient epochs
of launch times.

The above estimates are based on patching together three different point-mass,
two-body models (Moon + TP, Earth + TP, and finally Sun + TP). For the
purposes of these simple estimates, we have also ignored the effect of the
Moon’s rotation on the launch velocities of particles as well as the
eccentricity of the lunar orbit and of Earth’s heliocentric orbit. In detail,
the orbital evolution of escaping lunar ejecta particles that enter Earth-like
heliocentric orbits is subject to strongly chaotic dynamics and is exceedingly
sensitive to initial conditions, as is well known in the three-body problem,
hence the need for the numerical approach that follows below. These
theoretical estimates provide a guide for the initial conditions of the lunar
ejecta that are to be explored with numerical simulations and a guide for the
analysis of the results.

###  5.2 Numerical Model

We explore the dynamical fates of lunar ejecta in a similar vein as has been
done for satellite ejecta in the Saturnian system (aD10, ). Our dynamical
model includes the eight major planets from Mercury to Neptune and the Moon,
and we use the IAS15 integrator within REBOUND (hRdS14, ). The Direct
predefined module in REBOUND was used to detect collisions with the massive
bodies. An initial step size of 1.2 days was used and the step–size control
parameter was set to its default value ($\epsilon=10^{-9}$; this assures
machine precision for long time orbit propagations of $10^{10}$ orbital
periods (hRdS14, )). The length of the main set of simulations was 5,000
years; this is sufficiently long to explore the details of possible co-orbital
outcomes as a first study of the proof-of-concept for the lunar-ejecta
hypothesis for the origin of Kamo‘oalewa. In a second set of simulations, we
extended the simulation time up to 100,000 years for those particles
exhibiting Kamo‘oalewa-like dynamical behavior. Running these simulations for
much larger time spans is computationally prohibitive because, in order to
detect the QS and HS dynamical states and transitions between such co-orbital
configurations, it is necessary to have high cadence outputs ($\sim$ 1 output
per year); this places high demands on data storage.

The initial conditions for the planets are obtained from the JPL Horizons
system at epoch J2452996, i.e., 22 December 2003. In this initial exploration,
we adopt the Copernican principle (Peacock_1998, ; Gott93, ) that the current
time is not special, and is not unrepresentative of lunar ejecta launch
conditions at any random time in the geologically recent past. This assumption
can be tested in the future by exploring different initial epochs that sample
different initial conditions of the planets — especially of Earth’s
eccentricity — on secular timescales.

The initial conditions for the test particles are generated through three
parameters: the angle $\theta_{1}$ between the line joining the center of the
Moon to the launch site on the lunar surface and the line joining the center
of the Moon to center of Earth, the angle $\theta_{2}$ between the launch
velocity vector and the local normal vector at the launch site on the lunar
surface, and the speed of launch ${v}_{\text{L}}$. For simplicity, we consider
a two dimensional projection of the Moon’s surface onto the ecliptic,
illustrated in Fig. 2. Accordingly, the position of the TP is completely
specified by $\theta_{1}$ (and the Moon’s radius), while its initial velocity
(relative to the lunar surface) is determined by the magnitude and direction
of the specified relative velocity (${v}_{\text{L}}$ and $\theta_{2}$,
respectively); according to the projection made, this relative velocity has no
component perpendicular to the ecliptic plane. The angles $\theta_{1}$ and
$\theta_{2}$ range from $(0^{\circ},360^{\circ})$ and
$(-90^{\circ},90^{\circ})$, respectively. The values of ${v}_{\text{L}}$ were
chosen between 2.4 and 6.0 km s-1, with the lower bound corresponding to the
lunar-escape speed and the upper bound being the limit of ejection velocities
reported in numerical simulation studies of lunar cratering events (nAvS08, ).
It should also be noted that the frequency of ejection velocities decreases
rapidly as ${v}_{\text{L}}$ increases (nAvS08, ), discouraging the exploration
of larger values of ${v}_{\text{L}}$.

Following the guidance from §5.1, four launch sites were sampled on the Moon’s
surface. These four sites are representative of each of the four hemispheres
of the Moon: the near-side, far-side, leading side, and trailing side. The
locations of these are shown in Fig. 2. At each location, the launch speed was
varied systematically from 2.4 to 6.0 km s-1 (in increments of 0.1 km s-1),
and, for each speed, 100 particles were launched with different angles
$\theta_{2}$ (uniformly chosen along $-90^{\circ}$ and $90^{\circ}$). In
total, we launched 14,800 test particles.

In the simulations, we monitored for collisions with all the massive bodies.
In order to identify co-orbital outcomes, we visually examined the time series
of $a$ and $\Delta\lambda$. Rather than examining the time series of all
launched particles, this task is made easier by first projecting the evolution
in the $(a,e)$ plane and identifying those particles that appear in a rather
sparsely populated, narrow vertical zone in the semi-major axis range
$0.98-1.02$ au, as explained further in Section 2.

Acknowledgements We acknowledge an allocation of computer time from the UA
Research Computing High Performance Computing (HPC) at the University of
Arizona. The citations in this paper have made use of NASA’s Astrophysics Data
System Bibliographic Services. RM additionally acknowledges research support
from the program “Alien Earths” (supported by the National Aeronautics and
Space Administration under agreement No. 80NSSC21K0593) for NASA’s Nexus for
Exoplanet System Science (NExSS) research coordination network sponsored by
NASA’s Science Mission Directorate, and from the Marshall Foundation of
Tucson, AZ. AR acknowledges support by the Air Force Office of Scientific
Research (AFOSR) under Grant No. FA9550-21-1-0191.

Data availability Outcomes of the numerical simulations presented in this
paper are available at this publicly accessible permanent repository:
https://doi.org/10.5281/zenodo.8339513.

Code availability REBOUND can be freely downloaded from the developers’
webpage https://rebound.readthedocs.io. Codes and scripts used for the
numerical simulations may be requested to the corresponding author.

Competing interests The authors declare no competing interests.

Authors contributions JDCC wrote the codes and scripts used for the numerical
simulations, conducted the simulations, data processing, and data analysis. RM
composed the theoretical estimates and carried out related dynamical
calculations. RM and AJR designed the research project. All authors discussed
and critiqued interim and final results and contributed to writing the
manuscript.

## References

  *   * (1) Murray, C.D., Dermott, S.F.: Solar System Dynamics. Cambridge University Press, Cambridge (1999)
  * (2) Namouni, F., Christou, A.A., Murray, C.D.: Coorbital dynamics at large eccentricity and inclination. Phys. Rev. Lett. 83, 2506–2509 (1999)
  * (3) Gladman, P., Michel, P., Froeschlé, C.: The Near-Earth object population. Icarus 146, 176–189 (2000)
  * (4) Kaplan, M., Cengiz, S.: Horseshoe co-orbitals of Earth: current population and new candidates. Mon. Not. R. Astron. Soc. 496, 4420–4432 (2020)
  * (5) de la Fuente Marcos, C., de la Fuente Marcos, R.: Using Mars co-orbitals to estimate the importance of rotation-induced YORP break-up events in Earth co-orbital space. Mon. Not. R. Astron. Soc. 501, 6007–6025 (2021)
  * (6) Di Ruzza, S., Pousse, A., Alessi, E.M.: On the co-orbital asteroids in the solar system: medium-term timescale analysis of the quasi-coplanar objects. Icarus 390, 115330–20 (2023)
  * (7) de la Fuente Marcos, C., de la Fuente Marcos, R.: Asteroid 2014 $OL_{339}$: yet another Earth quasi-satellite. Mon. Not. R. Astron. Soc. 445, 2985–2994 (2014)
  * (8) Dermawan, B.: Temporal Earth coorbital types of asteroid 2016 HO3. J. Phys. Conf. Ser. 1523, 012019–4 (2019)
  * (9) Rezky, M., Soegiartini, E.: The orbital dynamics of asteroid 469219 Kamo‘oalewa. J. Phys. Conf. Ser. 1523, 012019–6 (2020)
  * (10) Fenucci, M., Novaković, B.: The role of the Yarkovsky effect in the long-term dynamics of asteroid (469219) Kamo‘oalewa. Astron. J. 162, 227–11 (2021)
  * (11) Sharkey, B.N.L., Reddy, V., Malhotra, R., et al.: Characterizing Earth quasi-satellite (469219) 2016 HO3 Kamo‘oalewa. Nature Commun. Earth Environ. 2, 231–7 (2021)
  * (12) de la Fuente Marcos, C., de la Fuente Marcos, R.: Asteroid (469219) 2016 HO3, the smallest and closest Earth quasi-satellite. Mon. Not. R. Astron. Soc. 462, 3441–3456 (2016)
  * (13) Brasser, R., Innanen, K.A., Connors, M., et al.: Transient co-orbital asteroids. Icarus 171, 102–109 (2004)
  * (14) Wajer, P.: 2002 AA29: Earth’s recurrent quasi-satellite. Icarus 200, 147–153 (2009)
  * (15) Venigalla, C., Baresi, N., Aziz, J.D., et al.: Near-Earth Asteroid Characterization and Observation (NEACO) mission to asteroid (469219) 2016 HO3. J. Spacecr. Rockets 56, 1121–1136 (2019)
  * (16) Jin, W.T., Li, F., Yan, J.G., et al.: A simulated global GM estimate of the asteroid 469219 Kamo‘oalewa for China’s future asteroid mission. Mon. Not. R. Astron. Soc. 493, 4012–4021 (2020)
  * (17) de la Fuente Marcos, C., de la Fuente Marcos, R.: (309239) 2007 RW10: a large temporary quasi-satellite of Neptune. Astron. Astrophys. 545, 9–4 (2012)
  * (18) Morais, M.H.M., Morbidelli, A.: The population of near-Earth asteroids in coorbital motion with the Earth. Icarus 160, 1–9 (2002)
  * (19) Gladman, B.J., Burns, J.A., Duncan, M.J., Levison, H.F.: The dynamical evolution of lunar impact ejecta. Icarus 118, 302–321 (1995)
  * (20) Granvik, M., Morbidelli, A., Jedicke, R., Bolin, B., Bottke, W.F., Beshore, E., Vokrouhlický, D., Nesvorný, D., Michel, P.: Debiased orbit and absolute-magnitude distributions for near-Earth objects. Icarus 312, 181–207 (2018) arXiv:1804.10265 [astro-ph.EP]. https://doi.org/10.1016/j.icarus.2018.04.018
  * (21) Nesvorný, D., Deienno, R., Bottke, W.F., Jedicke, R., Naidu, S., Chesley, S.R., Chodas, P.W., Granvik, M., Vokrouhlický, D., Brož, M., Morbidelli, A., Christensen, E., Shelly, F.C., Bolin, B.T.: NEOMOD: A New Orbital Distribution Model for Near-Earth Objects. Astron. J. 166(2), 55 (2023) arXiv:2306.09521 [astro-ph.EP]. https://doi.org/10.3847/1538-3881/ace040
  * (22) Bottke, W.F. Jr., Nolan, M.C., Melosh, H.J., et al.: Origin of the Spacewatch small Earth-approaching asteroids. Icarus 122, 406–427 (1996)
  * (23) Small Body Database Lookup, 2020 PN1. https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=54050997 Accessed 2023-06-30
  * (24) Small Body Database Lookup, 2020 PP1. https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=2020%20PP1 Accessed 2023-06-30
  * (25) Small Body Database Lookup, 2020 KZ2. https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=2020%20KZ2 Accessed 2023-06-30
  * (26) Ito, T., Malhotra, R.: Asymmetric impacts of near-Earth asteroids on the Moon. Astron. Astrophys. 519, 63–9 (2010)
  * (27) Gallant, J., Gladman, B., Ćuk, M.: Current bombardment of the Earth-Moon system: Emphasis on cratering asymmetries. Icarus 202, 371–382 (2009)
  * (28) Liu, T., Luther, R., Manske, L., Wünnemann, K.: Melt production and ejection at lunar intermediate-sized impact craters: where is the molten material deposited? In: European Planetary Science Congress, pp. 2022–402 (2022)
  * (29) Melosh, H.J.: High-velocity solid ejecta fragments from hypervelocity impacts. International Journal of Impact Engineering 5(1), 483–492 (1987). Hypervelocity Impact Proceedings of the 1986 Symposium
  * (30) Singer, K.N., Jolliff, B.L., McKinnon, W.B.: Lunar secondary craters and estimated ejecta block sizes reveal a scale-dependent fragmentation trend. J. Geophys. Res. Planets 125, 2019–00631327 (2020)
  * (31) Melosh, H.J.: An Empirical Function Linking Impact Ejecta Fragment Size and Velocity. In: 51st Annual Lunar and Planetary Science Conference. Lunar and Planetary Science Conference, p. 2587 (2020)
  * (32) Stöffler, D., Ryder, G., Ivanov, B.A., Artemieva, N.A., Cintala, M.J., Grieve, R.A.: Cratering history and lunar chronology. Reviews in Mineralogy and Geochemistry 60(1), 519–596 (2006)
  * (33) Near-Earth Object Surveyor Mission (2022). https://neos.arizona.edu/ Accessed 2023-06-30
  * (34) Shen, Y., Tremaine, S.: Stability analysis of Earth co-orbital objects. Astron. J. 136, 2453–2467 (2008)
  * (35) Pousse, A., Robutel, P., Vienne, A.: On the co-orbital motion in the planar restricted three-body problem: the quasi-satellite motion revisited. Celest. Mech. Dyn. Astr. 128, 383–407 (2017)
  * (36) Voyatzis, G., Antoniadou, K.I.: On quasi-satellite periodic motion in asteroid and planetary dynamics. Celest. Mech. Dyn. Astr. 130, 59–18 (2018)
  * (37) Hénon, M.: Numerical exploration of the restricted problem. V. Hill’s case: periodic orbits and their stability. Astron. Astrophys. 1, 223–238 (1969)
  * (38) Scott, C.J., Spencer, D.B.: Calculating transfer families to periodic distant retrograde orbits using differential correction. J. Guid. Control. Dyn. 33, 1592–1605 (2010)
  * (39) Lara, M.: Nonlinear librations of distant retrograde orbits: a perturbative approach—the Hill problem case. Nonlinear Dyn. 93, 2019–2038 (2018)
  * (40) Mikkola, S., Innanen, K., Wiegert, P., et al.: Stability limits for the quasi-satellite orbit. Mon. Not. R. Astron. Soc. 369, 15–24 (2006)
  * (41) Sidorenko, V.V., Neishtadt, A.I., Artemyev, A.V., Zelenyi, L.M.: Quasi-satellite orbits in the general context of dynamics in the 1:1 mean motion resonance: perturbative treatment. Celest. Mech. Dyn. Astr. 120, 131–162 (2014)
  * (42) Morais, M.H.M., Namouni, F.: Periodic orbits of the retrograde coorbital problem. Mon. Not. R. Astron. Soc. 490, 3799–3805 (2018)
  * (43) Oshima, K.: The roles of $L_{4}$ and $L_{5}$ axial orbits in transport among co-orbital orbits. Mon. Not. R. Astron. Soc. 480, 2945–2952 (2018)
  * (44) Qi, Y., Qiao, D.: Stability analysis of Earth co-orbital objects. Astron. J. 163, 211–14 (2022)
  * (45) Wajer, P.: Dynamical evolution of Earth’s quasi- satellites: 2004 GU$9$ and 2006 FV$35$. Icarus 209, 488–493 (2010)
  * (46) Christou, A.A., Asher, D.J.: A long-lived horseshoe companion to the Earth. Mon. Not. R. Astron. Soc. 414, 2965–2969 (2011)
  * (47) Connors, M.: A Kozai-resonating Earth quasi-satellite. Mon. Not. R. Astron. Soc. 437, 85–89 (2014)
  * (48) Christou, A.A., Georgakarakos, N.: Long-term dynamical survival of deep Earth co-orbitals. Mon. Not. R. Astron. Soc. 507, 1640–1650 (2021)
  * (49) Barrabé, E., Ollè, M.: Invariant manifolds of $L_{3}$ and horseshoe motion in the restricted three-body problem. Nonlinearity 19, 2065–2089 (2008)
  * (50) Qi, Y., Qiao, D.: Co-orbital transition of 2016 HO3. Astrodyn. ( in press, 2022)
  * (51) Todorović, N., Wu, D., Rosengren, A.J.: The arches of chaos in the Solar System. Sci. Adv. 6, 1313–6 (2020)
  * (52) Wieczorek, M.A., Le Feuvre, M.: Did a large impact reorient the moon? Icarus 200, 358–366 (2009)
  * (53) Aharonson, O., Goldreich, P., Sari, R.: Why do we see the man in the Moon? Icarus 219, 241–243 (2012)
  * (54) Dobrovolskis, A.R., Alvarellos, J.L., Zahnle, K.J., Lissauer, J.J.: Exchange of ejecta between Telesto and Calypso: Tadpoles, horseshoes, and passing orbits. Icarus 210, 436–445 (2010)
  * (55) Rein, H., Spiegel, D.S.: IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. Mon. Not. R. Astron. Soc. 446, 1424–1437 (2014)
  * (56) Peacock, J.A.: Cosmological Physics, (1999)
  * (57) Gott, I. J. R.: Implications of the Copernican principle for our future prospects. Nature 363(6427), 315–319 (1993)
  * (58) Artemieva, N.A., Shuvalov, V.V.: Numerical simulation of high-velocity impact ejecta following falls of comets and asteroids onto the Moon. Sol. Syst. Res. 42, 329–334 (2008)