A Collision Mechanism for the Removal of Earth's Trojan Asteroids

The restricted circular three-body problem predicts five points of equilibrium called the Lagrange points. Three of these points (L1, L2, and L3) are colinear with the primary and secondary masses, and the remaining two (L4 and L5) lead and trail the secondary mass by 60°, respectively. The three colinear points are all saddle points of the potential, and are therefore unstable. However, L4 and L5 are points of linearly stable equilibrium (as long as the secondary/primary mass ratio μ < 0.039; see Murray & Dermott 2000), and, as such, can host long-term stable populations, colloquially known as Trojan asteroids.

Because Trojan asteroids are in the 1:1 mean-motion resonance with their host planet, they can remain stable for billions of years. As a result, they have the potential to be of primordial origin and can provide important probes of the solar system's dynamical and chemical history. The Jupiter Trojans are of considerable interest, as their quantity and chemical composition may encode sensitive information about the composition of the protoplanetary disk and the migration history of the giant planets. Understanding this population is of such high importance that NASA's Lucy spacecraft will fly by seven Jupiter Trojans over the next decade, four of which are in binary systems (Levison et al. 2021). Neptune Trojans encode similar information and are now being studied in detail, with characterization efforts including a Cycle 1 James Webb Space Telescope program to obtain near-infrared reflectance spectra (Markwardt et al. 2021).

Only two Earth Trojans (ETs) have been found to date, and both are transient; simulations show that they are stable over timescales much shorter than the age of the solar system (Connors et al. 2011; Hui et al. 2021). The prospect of discovering long-term stable ETs is enticing, as their proximity and low delta-v with respect to Earth would make a spacecraft mission imminently possible. Additionally, it would likely be possible to return samples of ETs. This feature would be extremely valuable for studying the composition of the protoplanetary disk if the objects are indeed of primordial origin. If such Trojans are not of primordial origin, the samples would still be useful for gleaning insight about the reservoirs that populate near-Earth orbits. In light of the potential scientific payoff of discovering ETs, several surveys have been dedicated to searching for such objects (Whiteley & Tholen 1998; Cambioni et al. 2018; Yoshikawa et al. 2018; Markwardt et al. 2020; Lifset et al. 2021). While these searches have failed to discover any long-term stable ETs, they have placed significant observational upper limits on the possible population. Specifically, the number of ETs with sizes larger than approximately 400 meters is limited to at most 100 objects (Markwardt et al. 2020; Lifset et al. 2021).

These observations are seemingly in tension with numerical simulations showing that ETs can be stable over billion-year timescales (see, e.g., Tabachnik & Evans 2000; Ćuk et al. 2012; Marzari & Scholl 2013; Malhotra 2019; Zhou et al. 2019; Christou & Georgakarakos 2021). Furthermore, while it has been shown that the Yarkovsky effect (Bottke et al. 2006; Marzari & Scholl 2013; Zhou et al. 2019) can remove the small ETs (smaller than tens of kilometers) from resonance, the effect is not strong enough to remove the larger bodies on timescales comparable to the age of the solar system. Current theories therefore suggest that long-term stable ETs can exist. To date, any such object has eluded discovery.

Although the existence of any large ETs has not been ruled out by observational surveys, the continued absence of even a single long-term stable ET merits the consideration of a gap in our understanding of their dynamics. In this work, we consider how collisions between asteroids and the Earth, such as those hypothesized to have occurred during the Moon-forming impact (see, e.g., Canup 2012; Ćuk & Stewart 2012) or while Earth was acquiring its late veneer (see, e.g., Kleine 2011; Raymond et al. 2013; Brasser et al. 2016), would affect an initially stable ET population. In Section 2 we derive the circumstances of the collisions between the Earth and the asteroids that may have been responsible for delivering its late veneer. In Section 3, we describe our numerical experiments in which we simulated the bombardment of Earth by large asteroids. In Section 4 we present the results of our simulations and analyze their outcomes. In Section 5, we summarize our results and interpret them in the broad context of the possible ET population. For completeness, Appendix A provides a quantitative assessment of some of our approximations and Appendix B presents an order-of-magnitude estimate for the orbital changes required to remove Trojans from bound states. Finally, we consider perturbations of Earth's orbit due to close encounters that do not result in collisions (Appendix C) and direct collisions between incoming asteroids and Trojans (Appendix D).

In this section we compute the distributions of relative velocity and orientation for potential Earth impactors (i.e., asteroids on Earth-crossing orbits). We use the resulting distributions to compute the change in Earth's momentum for our simulations of asteroid–Earth collisions in Section 3.

2.1. Scales and Ordering

Consider an asteroid on an Earth-crossing orbit. Here the asteroid has mass m, the Earth has mass M_⊕, and the Sun has mass M_⊙. The masses obey the ordering

$$\begin{eqnarray}&&m\ll {M}_{\oplus }\ll {M}_{\odot },\end{eqnarray} \tag{ 1 }$$

so that the orbits can be considered in the hierarchical limit. Using the mass of Earth, the mass ratio μ = M_⊕/M_⊙ ∼3 × 10⁻⁶.

As a first approximation, we take Earth to have a circular orbit with semimajor axis a_⊕ = 1 au. The potential impactor has initial orbital elements (a, e) and is considered to have the same orbital plane as Earth (i = 0). The Hill sphere of Earth is given by

$$\begin{eqnarray}&&{R}_{H}={a}_{\oplus }{\left(\displaystyle \frac{{M}_{\oplus }}{3{M}_{\odot }}\right)}^{1/3}\approx 0.01\ \mathrm{au},\end{eqnarray} \tag{ 2 }$$

whereas the sphere of influence has the form

$$\begin{eqnarray}&&{R}_{S}={a}_{\oplus }{\left(\displaystyle \frac{{M}_{\oplus }}{{M}_{\odot }}\right)}^{2/5}\approx 0.006\ \mathrm{au}.\end{eqnarray} \tag{ 3 }$$

Since the orbital length scales obey the ordering

$$\begin{eqnarray*}&&{R}_{S}\sim {R}_{H}\ll {a}_{\oplus }\lt a,\end{eqnarray*}$$

we can use intermediate asymptotics to describe the dynamics of the bodies that collide with Earth. In other words, the orbit of the asteroid has a Keplerian orbit around the Sun until it reaches the vicinity of Earth. Once the asteroid enters the Earth's sphere of influence, its subsequent dynamics are governed by the gravitational potential of Earth. We can thus find the velocity vector of the asteroid at the point where it crosses the orbit of Earth. This velocity vector then provides the outer boundary condition for the inner problem where the asteroid collides or interacts with Earth.

2.2. Relative Velocity for Earth Encounters

We now calculate the speed of the asteroid with respect to Earth at the point of orbit crossing. The asteroid orbit has specific energy and angular momentum, given by

$$\begin{eqnarray}&&E=-\displaystyle \frac{{{GM}}_{\odot }}{2a}\qquad \mathrm{and}\qquad {J}^{2}={{GM}}_{\odot }a(1-{e}^{2}).\end{eqnarray} \tag{ 4 }$$

The azimuthal velocity component of the asteroid orbit is determined by conservation of angular momentum so that

$$\begin{eqnarray}&&{v}_{\phi \odot }^{2}=\displaystyle \frac{{J}^{2}}{{r}^{2}}=\displaystyle \frac{{{GM}}_{\odot }a}{{r}^{2}}(1-{e}^{2}).\end{eqnarray} \tag{ 5 }$$

Note that the azimuthal velocity is defined here in the solar reference frame, as indicated by the subscript. The radial velocity of the asteroid with respect to the Sun is then given by

$$\begin{eqnarray}\begin{array}{rcl}{v}_{r}^{2} & = & 2E-\displaystyle \frac{{J}^{2}}{{r}^{2}}+\displaystyle \frac{2{{GM}}_{\odot }}{r}=\displaystyle \frac{2{{GM}}_{\odot }}{r}\\ & & -\displaystyle \frac{{{GM}}_{\odot }a}{{r}^{2}}(1-{e}^{2})-\displaystyle \frac{{{GM}}_{\odot }}{a}.\end{array}\end{eqnarray} \tag{ 6 }$$

Because Earth is assumed to have a circular orbit, it has no radial velocity component, so that the radial velocity of the asteroid is the same in both reference frames. In this approximation, the speed of the asteroid with respect to Earth is thus given by

$$\begin{eqnarray}&&{v}_{\infty }^{2}={v}_{r}^{2}+{\left({v}_{\phi \odot }-{v}_{\oplus }\right)}^{2}={v}_{r}^{2}+{v}_{\phi \odot }^{2}+{v}_{\oplus }^{2}-2{v}_{\phi \odot }{v}_{\oplus },\end{eqnarray} \tag{ 7 }$$

where ${v}_{\oplus }={({{GM}}_{\odot }/{a}_{\oplus })}^{1/2}\approx 30$ km s⁻¹ is the orbital speed of Earth, and where all quantities are evaluated at r = a_⊕. ³ Using the results obtained above, we find

$$\begin{eqnarray}&&{v}_{\infty }^{2}=\displaystyle \frac{{{GM}}_{\odot }}{{a}_{\oplus }}\left[3-\displaystyle \frac{{a}_{\oplus }}{a}-2{\left(\displaystyle \frac{a}{{a}_{\oplus }}\right)}^{1/2}{(1-{e}^{2})}^{1/2}\right].\end{eqnarray} \tag{ 8 }$$

We see that the orbital speed of Earth, v_⊕, sets the velocity scale for this collision problem. Let us define the dimensionless variables

$$\begin{eqnarray}&&{u}_{\infty }\equiv \displaystyle \frac{{v}_{\infty }}{{v}_{\oplus }}\qquad \mathrm{and}\qquad \xi \equiv \displaystyle \frac{a}{{a}_{\oplus }}\end{eqnarray} \tag{ 9 }$$

so that the expression for the hyperbolic speed at infinity reduces to the form

$$\begin{eqnarray}&&{u}_{\infty }^{2}\,=\,3-\displaystyle \frac{1}{\xi }-2{\xi }^{1/2}{(1-{e}^{2})}^{1/2}.\end{eqnarray} \tag{ 10 }$$

The largest possible speed is thus ${u}_{\infty }=\sqrt{3}$ or ${v}_{\infty }={v}_{\oplus }\sqrt{3}\approx 52$ km s⁻¹. ⁴ This case corresponds to a distant asteroid orbit (large a) with small perihelion (eccentricity e close to unity) so that the asteroid velocity is nearly radial when it encounters Earth. For more typical cases, we expect perihelion q = ξ(1 − e) ∼ 1 in dimensionless units, and ξ ∼ 2, so that u_∞ ∼ 0.2 and v_∞ ∼ 6 km s⁻¹. Note that the expressions derived above also apply for orbits with a < a_⊕ (although such objects do not originate from the asteroid belt).

2.3. Distributions of Impact Variables

Using the results derived above, we can find the distribution of relative speeds by assuming a simple model for the initial population of asteroid impactors. We can also determine the distribution of angles at which the asteroids intercept the orbit of Earth. A great deal of previous work has quantified the distributions of speeds (and other properties) for near-Earth asteroids for purposes of planetary defense (e.g., see Steel 1998; Bottke et al. 2000; Jeffers et al. 2001; and many others). These distributions are consistent with the model developed below (although the present-day and primordial distributions could be different).

Let us consider the case where the asteroids have a surface density distribution of the form σ ∝ 1/a, so that the distribution of the semimajor axes of the asteroids is uniformly distributed within some range, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{da}}=\mathrm{constant}\qquad \mathrm{for}\qquad {a}_{1}\leqslant a\leqslant {a}_{2}.\end{eqnarray} \tag{ 11 }$$

For the sake of definiteness, here we assume a₁ = 1 au and a₂ = 5 au, so we are sampling the full annulus between the current orbits of Earth and Jupiter. For a given value of a, the eccentricity must be larger than a minimum value given by

$$\begin{eqnarray}&&{e}_{\min }=1-\displaystyle \frac{{a}_{\oplus }}{a}.\end{eqnarray} \tag{ 12 }$$

If this condition is not met, then the asteroid can not cross the orbit of Earth. The eccentricity is unlikely to be much larger than that corresponding to an Earth-crossing orbit: it takes interactions to pump up the eccentricity, and continued orbit crossings will eventually lead to large orbital perturbations. For the sake of definiteness, let us assume that the maximum eccentricity corresponds to perihelion of a_⊕/10, so that

$$\begin{eqnarray}&&{e}_{\max }=1-\displaystyle \frac{{a}_{\oplus }}{10a}.\end{eqnarray} \tag{ 13 }$$

We then assume that the distribution of eccentricities is uniform over the interval $[{e}_{\min },{e}_{\max }]$. With these choices for the distributions of the asteroid orbital elements (a, e), the resulting distribution of relative speeds, v_∞, is shown in Figure 1.

Zoom In Zoom Out Reset image size

We can also find the distribution of directions. Let θ be the angle of the incoming velocity (in the planetary reference frame) with respect to the velocity of Earth. In other words,

$$\begin{eqnarray}&&\cos \theta =\displaystyle \frac{{\boldsymbol{v}}\cdot {{\boldsymbol{v}}}_{\oplus }}{| {\boldsymbol{v}}| | {{\boldsymbol{v}}}_{\oplus }| }=\displaystyle \frac{({v}_{\phi \odot }-{v}_{\oplus })}{{v}_{\infty }}.\end{eqnarray} \tag{ 14 }$$

The resulting distributions of both relative speed and angular direction are shown in Figure 1. These distributions are used as the starting point for the numerical simulations discussed in the following section. Note that the distribution of angles has two peaks, corresponding to bodies impacting Earth on incoming and outgoing trajectories. Moreover, the two parts of distribution, centered on θ = π/2 and 3π/2, each include both positive and negative values of $\cos \theta$. In Appendix A we investigate the validity of the approximations made in this section.

The primary contribution of this work is an ensemble of numerical experiments where we tested the stability of a synthetic population of long-term stable ETs while Earth was stochastically bombarded by asteroids. For all of our numerical experiments, we used Rebound's WHFAST integrator (Rein & Tamayo 2015) and a model of the solar system that included the Sun, Venus, Earth, Mars, Jupiter, Saturn, Uranus, and Neptune. ⁵

We began by generating a synthetic population of 3021 ETs, each of which was dynamically stable for at least 10 Myr. This population represents the stable subset of an initial set of candidate ETs (which were integrated numerically to verify stability). The surviving objects were split roughly evenly between L4 and L5. We show the resulting population projected into the ecliptic plane in Figure 2. The dynamical properties of our synthetic population match those reported by Zhou et al. (2019). Note that due to the presence of other planets, these stable clouds are wider than one would expect by considering only the circular restricted three-body problem (e.g., Murray & Dermott 2000).

Zoom In Zoom Out Reset image size

Next we ran simulations to test the stability of our synthetic ETs as Earth underwent collisions with asteroids. We began each simulation by setting Earth's mass to 0.99 M_⊕, and integrating for a warm-up period of 100 yr with no collisions. After the warm-up period, we perturbed Earth stochastically by modifying its velocity as it would change due to a collision with an asteroid. ⁶ We randomly drew the impactor velocity and impact angle from the distributions shown in Figure 1. We chose the impactor mass by drawing a radius, R, from a distribution of the form N(R) ∝ R⁻² with values ranging from 500 to 2000 km, and assuming a mass density ρ = 3 g cm⁻³. While there are certainly many asteroids smaller than 500 km, they will deliver a negligible amount of mass to Earth compared to the mass delivered by a few large objects. Once we had perturbed Earth, we integrated for 10,000 more years (we found that this was enough time for the system to roughly settle back into equilibrium), tracking whether each test particle remained a Trojan, migrated onto a horseshoe orbit, or was ejected from a coorbital state altogether. If an object was no longer a coorbital with Earth, we removed it from the simulation. We continued this procedure until Earth's mass reached 1 M_⊕. ⁷ We repeated this process a total of 5000 times to build up our statistics. The numerical values used in these simulations were all chosen to be roughly consistent with the circumstances of Earth's late veneer impacts reported in the literature (see, e.g., Raymond et al. 2013; Brasser et al. 2016, and the references therein).

In this section we analyze the data from our simulations. We begin by examining the probabilities with which initially-Trojan objects remained on true Trojan orbits, or were perturbed onto horseshoe orbits. In Figure 3, we show box-and-whisker plots of the fraction of the initial population of 3021 ETs that ended the simulation in a Trojan or a horseshoe orbit, as a function of the number of impactors. Note that there are data points that exceed the vertical range we have chosen for Figure 3, but such events are relatively rare. After running many simulations an apparent trend emerges: for any number of impacts, the median surviving fraction of Trojans was less than a few percent. In addition, for a given number of impactors, the median number of objects surviving on horseshoe orbits was always on the order of 1%–2%. These results indicate that large asteroid impacts with Earth are efficient at clearing the coorbital population.

Zoom In Zoom Out Reset image size

This series of collisions is a stochastic process in which different combinations of the number, velocities, and impact angles of the impactors can lead to vastly different outcomes. Indeed, there were some simulations in which an appreciable fraction of the objects remained on Trojan orbits, and as such require a more nuanced analysis. We find that all such simulations were instances where the number of impactors was small, and where all of them had low relative velocities. In these cases, the impactors simply did not impart enough momentum to Earth to remove the majority of the Trojan population. We demonstrate this trend in Figure 4 by plotting the surviving fractions of Trojan and horseshoe orbits as a function of the net linear momentum of the impactors. As the net linear momentum imparted by impactors increases, the probability of many objects remaining on Trojan orbits decreases. As the number of impactors increases, the probability of drawing all impactors with low velocities decreases. This decreases the probability of a small net linear momentum, and thus decreases the frequency and severity of outliers. Of course, there is a limit to this logic. In the case that Earth gradually acquired its late veneer from millions of small rocks, the net linear momentum would converge to zero, and the Trojans would remain largely unaffected. In a similar vain, the most common outcome (by far) with 1–5 impactors was that no objects remained in Trojan orbits at the end of the simulation. This outcome occurs because with so few impactors, a large net linear momentum is far more probable than a small net linear momentum (see the distributions in Figure 1).

Zoom In Zoom Out Reset image size

Finally, in Figure 5 we show the distribution of changes in Earth's semimajor axis and eccentricity during our simulations. We find that Earth's semimajor axis never changes by more than a few percent, and its eccentricity remains reasonably small—always within the range of predicted secular oscillations for the present-day solar system. The distribution of changes in semimajor axis is nearly symmetric about δa = 0. In contrast, the distribution of changes in eccentricity is slightly skewed toward negative values (decreases in eccentricity, i.e., increases in angular momentum). Note that we began our simulations with Earth's a = 1 au and e = 0.0167. The exact starting values are not important for this study, as long as they are in the vicinity of today's observed values. If the starting eccentricity is close to zero, however, then the changes in eccentricity (while still small) must be skewed toward the positive. Finally, we note that our choice of orbital elements for the incoming rocks includes specific angular momenta that are both larger and smaller than that of different choices of these distributions could lead to somewhat different changes in the orbit of Earth (see, e.g., Fernandez & Ip 1984).

Zoom In Zoom Out Reset image size

In this paper we have examined the stability of the ETs in the face of collisions between Earth and massive asteroids. The main result of this work is the demonstration that large impacts to Earth can readily remove Trojans from resonance. Using an ensemble of 5000 numerical simulations, we have demonstrated that if Earth acquired the final 1% of its present-day mass from a handful of asteroidal impactors, the population of ETs would be reduced by more than 99% (less than 1% remaining; see Figure 3). Due to the wide distribution of impactor parameters and the chaotic nature of the system, there are some instances in which a significant fraction of the original Trojans can remain on Trojan orbits. However, such instances are rather rare. In addition, we find that the population of objects on horseshoe orbits exceeds that of the Trojans, provided that the total number of collisions remains small. For the scenarios explored in this work, on average approximately 2% of the primordial Trojan population ends up on horseshoe orbits. ⁸ As long as the number of impactors remains small, the number of Trojans expected to survive is at least an order of magnitude smaller than the horseshoe objects. ⁹

We also analyzed our simulation data for correlations between the impactor parameters and the surviving fraction of Trojans. First we find that as the number of impactors increases to more than of order ∼10, the surviving fraction of Trojans increases. This finding makes sense intuitively: if we consider the limiting case in which Earth acquired the final 1% of its mass by having 10²⁸ grains of sand slowly and isotropically scattered onto its surface, the resulting impulse would be a negligible and the Trojans would remain undisturbed. In a similar vein, we find a trend between the maximum retention fraction of the Trojan population and the net linear momentum imparted to Earth by the collisions (see Figure 4). As the net linear momentum increases, the upper bound on the surviving fraction of Trojans decreases. We found that no Trojans survived in any simulation with net linear momentum Δp ≳ 0.25 M_⊕ km s⁻¹ (see also Appendix B). We also analyze the changes in Earth's semimajor axis and eccentricity at the end of the simulations to make sure that its orbit did not become unreasonably excited. The changes to the orbital elements are relatively small (see Figure 5), with Δa ∼ 0.01 au and Δe ∼ 0.005.

This set of results has important implications for our understanding of the fate of primordial ETs. While numerical simulations show that stable ETs can exist, dedicated searches still have not discovered any such objects. It has been noted that the Yarkovsky effect can remove small Trojans from resonance on the timescale of the age of the solar system, but it can not explain the apparent absence of large objects. In contrast, the collision mechanism explored in this paper is independent of the size of the Trojans. Moreover, this collision mechanism is of particular interest because previous work suggests that Earth has indeed experienced large impact events after its formation: one example being the impact that formed the Moon (e.g., Canup 2012), and another a series of impacts in which Earth acquired its late veneer (e.g., Raymond et al. 2013). The Moon is thought to have formed as the result of a Mars-sized object colliding with the primordial Earth. Our simulations show that such a collision is likely to have cleared out the primordial ETs (though perhaps not if the collision was relatively slow, thus imparting a rather small change in net linear momentum). However, a collision this large would probably have reinjected some material into Trojan resonances (note that similar ideas have been put forth in the context of Martian Trojans, e.g., Polishook et al. 2017). On the other hand, the collisions that gave Earth its late veneer are hypothesized to have occurred well after the Moon-forming collision. If the late veneer was comprised of a small number of large impactors, as favored by recent models (Raymond et al. 2013; Genda et al. 2017), any reinjected material is likely to have been cleared. Furthermore, the impacts comprising the late veneer were probably not large enough to repopulate L4 and L5 with ejecta (especially bodies large enough to survive the Yarkovsky effect), but this complication should be explored in future work. Additionally, because the late veneer is thought to have occurred late in the formation of the inner solar system, there was probably not a large population of ambient bodies to repopulate the Trojan resonances (Brasser et al. 2016). Despite the large parameter space in which dynamically stable ETs can exist, it is possible that all or most such objects have been removed due to asteroids impacting Earth.

We note that this work has limitations and can be improved in the future. Given a considerable amount of additional computing, future studies should by more thoroughly explore the distributions of impactor parameters. In particular, changing the size, mass, and eccentricity distributions of the impactors could influence the results. For example, a shallower size distribution would mean that fewer impacts are necessary to deliver the required mass to Earth, while a steeper size distribution would mean that more impacts are necessary. In this vein, the total mass of the impactors is also important, and should be considered in future studies. We also made the approximation that the collisions were coplanar with Earth's orbit, while in reality they would have some distribution of inclinations. One should also address the possibility of impactors coming from outside of the asteroid belt. In this work we considered only impactors on orbits with a ∈ [1, 5] au. Orbits with larger semimajor axes corresponding to the giant planet region or the Kuiper belt would have marginally larger relative velocities when entering Earth's Hill sphere. Such collisions would deliver larger impulses to Earth than their low-a counterparts, and would thus be more detrimental to the survival of the Trojan populations. However, this correction is rather small: a body with a → ∞ can achieve a maximum velocity of only ∼52 km s⁻¹. We also did not model the effect of close encounters between Earth and potential impactors, or impacts directly onto the Trojan asteroids. We make analytic estimates of the relative sizes of these effects in Appendices C and D, respectively. We find that the expected net velocity perturbation to Earth from close encounters is roughly equivalent to that expected from direct collisions. The perturbations from close encounters will be orthogonal to those from direct collisions, so including close encounters in the model should enhance the net perturbation to Earth by a factor of approximately $\sqrt{2}$. Direct collisions with Trojans, on the other hand, are a less likely source of orbital disruption than perturbations to Earth's orbit by an order of magnitude. These considerations imply that the present work likely represents a lower limit on the disturbance to the Trojan population. Of course, the most important additional work is observational: future observational surveys will either detect a (presumably small) population of ETs or place increasingly tight limits on their existence. With sufficiently tight constraints on the population of primordial ETs, future simulations can in turn constrain the properties of the primordial asteroid belt.

We thank Renu Malhotra for useful discussions, and two anonymous referees for useful comments. This material is based upon work supported by the National Aeronautics and Space Administration under grant No. NNX17AF21G issued through the SSO Planetary Astronomy Program and by the National Science Foundation under grant No. AST-2009096.

In this section, we assess the severity of the approximations used in Section 2. In particular, we have assumed that the orbit of the asteroid only depends on the solar gravitational potential, and we have assumed that the orbit of Earth is circular.

In this approximation, we assume that the asteroid follows an elliptical orbit around the Sun, and then encounters Earth at the point where the orbits intersect. In this phase, the distance, d_⊙, between the asteroid and Earth in the frame of the solar system thus vanishes (d_⊙ → 0) at the intersection. The approximation then switches over to a reference frame where the planet (Earth) is at rest, and the asteroid is assumed to have speed v_∞ and distance d_⊙ → ∞ . We are thus using the method of matched intermediate asymptotics, where the inner boundary condition for the outer problem (d_⊙ → 0) represents the outer boundary condition for the inner problem (d_⊙ → ∞). In order to test the consistency of this approach, we need to determine how much the gravitational potential of Earth affects the incoming speed, v_∞, calculated previously.

The boundary between the inner problem and the outer problem is delineated by the Hill sphere, as given by Equation (2). At this Hill sphere boundary, the additional velocity that the asteroid would gain from the gravitational potential well of Earth is given approximately by

$$\begin{eqnarray}\begin{array}{rcl}{v}_{x}^{2} & = & \displaystyle \frac{2{{GM}}_{\oplus }}{{R}_{H}}=\displaystyle \frac{2{{GM}}_{\oplus }}{{a}_{\oplus }}{\left(\displaystyle \frac{3{M}_{\odot }}{{M}_{\oplus }}\right)}^{1/3}\\ & = & {3}^{1/3}\displaystyle \frac{2{{GM}}_{\odot }}{{a}_{\oplus }}{\left(\displaystyle \frac{{M}_{\oplus }}{{M}_{\odot }}\right)}^{2/3}.\end{array}\end{eqnarray} \tag{ A1 }$$

As shown previously, the incoming speed v_∞ is comparable to the orbital speed of Earth. If we let v_T denote the true speed (corrected for Earths gravity), then

$$\begin{eqnarray}&&{v}_{T}^{2}={v}_{\infty }^{2}+{v}_{x}^{2},\end{eqnarray} \tag{ A2 }$$

which in turn implies

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{T}}{{v}_{\infty }}={\left[1+\displaystyle \frac{{v}_{x}^{2}}{{v}_{\infty }^{2}}\right]}^{1/2}\approx 1+\displaystyle \frac{{v}_{x}^{2}}{2{v}_{\infty }^{2}}.\end{eqnarray} \tag{ A3 }$$

The size of the correction is thus given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}v}{v}=\displaystyle \frac{{v}_{T}-{v}_{\infty }}{{v}_{\infty }}\approx {3}^{1/3}{\left(\displaystyle \frac{{M}_{\oplus }}{{M}_{\odot }}\right)}^{2/3}\sim 0.0003.\end{eqnarray} \tag{ A4 }$$

If we use the sphere of influence from Equation (3), instead of the Hill radius, the above equation has the form

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}v}{v}=\displaystyle \frac{{v}_{T}-{v}_{\infty }}{{v}_{\infty }}\approx {\left(\displaystyle \frac{{M}_{\oplus }}{{M}_{\odot }}\right)}^{3/5}\sim 0.0005.\end{eqnarray} \tag{ A5 }$$

The correction due to the gravity of Earth is thus relatively small.

Another correction arises from our assumption that Earth has a circular orbit. In this case, the planet has a constant speed ${v}_{\oplus }^{2}={{GM}}_{\odot }/{a}_{\oplus }$. For small eccentricity, the variation in speed over the course of the orbit is given approximately by the expression

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\rm{\Delta }}v}{v}\right|}_{\max }={\left[\displaystyle \frac{1+e}{1-e}\right]}^{1/2}-1\approx e\approx 0.0167,\end{eqnarray} \tag{ A6 }$$

where we have used the observed eccentricity of Earth. In summary, both corrections are small, less than about 1%. The correction for the eccentricity of Earth's orbit is larger than the correction due to the Earths gravity.

Next we consider the crossing time: the relative asteroid speeds are of order v_∞ = 10 km s⁻¹ and the size of the Hill sphere is of order R_H = 10¹¹ cm. The crossing time of the asteroid is thus of order

$$\begin{eqnarray}&&{t}_{\mathrm{cross}}\approx \displaystyle \frac{{R}_{H}}{{v}_{\infty }}\sim {10}^{5}\ {\rm{s}}\,\sim \,1\ \mathrm{day}.\end{eqnarray} \tag{ A7 }$$

Over the crossing time (∼1 yr), the planet rotates through about one-three-hundredth of its orbit, so the reference frame of the planet does not rotate appreciably over the course of the encounter. We thus expect the effects of the rotating reference frame to be relatively small, but they can be quantified, as shown below. The encounter takes place in a rotating frame of reference, with rotation rate, Ω, given by the mean motion of Earth. In a rotating reference frame, the acceleration is given by

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{R}={{\boldsymbol{a}}}_{g}-2{\rm{\Omega }}\hat{z}\times \dot{{\boldsymbol{r}}}-{{\rm{\Omega }}}^{2}\hat{z}\times (\hat{z}\times {\boldsymbol{r}}),\end{eqnarray} \tag{ A8 }$$

where we have taken the angular velocity vector of the rotating frame to point in the $\hat{z}$ direction. On the right-hand side of the equation, the first term is the acceleration due to the gravity of Earth, the second term is the Coriolis acceleration, and the third term is the centrifugal acceleration. In order of magnitude, we can find the relative size of the Coriolis term and the centrifugal term,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{a}_{\mathrm{cor}}}{{a}_{g}} & = & 2{\left(\displaystyle \frac{{M}_{\odot }{r}^{3}}{{M}_{\oplus }{a}_{\oplus }^{3}}\right)}^{1/2}=\displaystyle \frac{2}{\sqrt{3}}{\left(\displaystyle \frac{r}{{R}_{H}}\right)}^{3/2}\quad \mathrm{and}\\ \displaystyle \frac{{a}_{\mathrm{cen}}}{{a}_{g}} & = & \left(\displaystyle \frac{{M}_{\odot }{r}^{3}}{{M}_{\oplus }{a}_{\oplus }^{3}}\right)=\displaystyle \frac{1}{3}{\left(\displaystyle \frac{r}{{R}_{H}}\right)}^{3}.\end{array}\end{eqnarray} \tag{ A9 }$$

At the Hill radius, both ratios are thus of order unity, but they both decrease rapidly as the asteroid penetrates the sphere and orbits closer to the planet. The centrifugal term is always smaller than the Coriolis term. We can thus estimate the size of this correction by using only the Coriolis term, and by noting that the change in velocity produced by this additional acceleration is given approximately by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}v & \sim & \displaystyle \int {a}_{\mathrm{cor}}{dt}\sim \displaystyle \int 2{\rm{\Omega }}{vdt}\sim 2{\rm{\Omega }}{dr}\lt 2{\rm{\Omega }}{R}_{H}\\ & = & 2{\left(\displaystyle \frac{{{GM}}_{\odot }}{{a}_{\oplus }}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{\oplus }}{3{M}_{\odot }}\right)}^{1/3}.\end{array}\end{eqnarray} \tag{ A10 }$$

We thus obtain the approximate bound

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}v}{v}\lesssim 2{\left(\displaystyle \frac{{M}_{\oplus }}{3{M}_{\odot }}\right)}^{1/3}\approx 0.03.\end{eqnarray} \tag{ A11 }$$

If we use the sphere of influence instead of the Hill radius, the above expression becomes

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}v}{v}\lesssim 2{\left(\displaystyle \frac{{M}_{\oplus }}{{M}_{\odot }}\right)}^{2/5}\approx 0.01.\end{eqnarray} \tag{ A12 }$$

The correction for the rotating frame of reference is thus small, but it is somewhat larger than the corrections outlined earlier due to the influence of Earth in the solar reference frame and the noncircular nature of its orbit.

Here we make some rough estimates of the energetic requirements for the removal of ETs from bound orbits. In this context, a bound orbit is one such that the body is bound to either the L4 or L5 point, and hence executes a Trojan orbit. Unbound orbits are no longer Trojans, but still remain bound to the Sun.

We work in the approximation of the circular restricted three-body problem, following the notation of Murray & Dermott (2000), and work in their dimensionless units where the mean motion of Earth is unity and the separation of the two majors bodies is also unity. In these units, the Jacobi constant has the form

$$\begin{eqnarray}&&{C}_{J}=({x}^{2}+{y}^{2})+2\left(\displaystyle \frac{1-\mu }{{r}_{1}}-\displaystyle \frac{\mu }{{r}_{2}}\right)-{v}^{2}.\end{eqnarray} \tag{ B1 }$$

The values of the Jacobi constants for the five Lagrangian equilibrium points are as follows:

$$\begin{eqnarray}&&{C}_{1}=3+{3}^{4/3}{\mu }^{2/3}-10\mu /3,\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&{C}_{2}=3+{3}^{4/3}{\mu }^{2/3}-14\mu /3,\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&{C}_{3}=3+\mu ,\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&{C}_{4}={C}_{5}=3-\mu .\end{eqnarray} \tag{ B5 }$$

We can thus define the difference in Jacobi constants between Lagrangian points:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{3-45}={C}_{3}-{C}_{4}\,=\,2\mu \end{eqnarray} \tag{ B6 }$$

and

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{12-45}=\displaystyle \frac{1}{2}\left({C}_{1}+{C}_{2}\right)-{C}_{4}={\left(9\mu \right)}^{2/3}-3\mu .\end{eqnarray} \tag{ B7 }$$

As a rough approximation, the velocity increment required to move an object from a Trojan orbit to a horseshoe orbit is then given by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}v & \sim & {\left({{\rm{\Delta }}}_{3-45}\right)}^{1/2}{v}_{\oplus }={\left(2\mu \right)}^{1/2}{v}_{\oplus }\\ & \approx & 0.002\,4{v}_{\oplus }\qquad \Rightarrow \qquad {\rm{\Delta }}v\sim 0.074\,\mathrm{km}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ B8 }$$

Similarly, the velocity increment required to make a Trojan become unbound from any libration is given approximately by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}v & \sim & {\left({{\rm{\Delta }}}_{12-45}\right)}^{1/2}{v}_{\oplus }={\left(9\mu \right)}^{1/3}{v}_{\oplus }\\ & \approx & 0.03{v}_{\oplus }\qquad \Rightarrow \qquad {\rm{\Delta }}v\sim 0.9\,\mathrm{km}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ B9 }$$

The above values represent upper bounds on the velocity increments required to dislodge the Trojans. Note that when an impact takes place on Earth, the immediate collision changes the velocity of Earth, so that the velocity of the Trojan in the Earth-centric frame changes accordingly. However, there is another effect: the collision also induces an eccentricity in Earth's orbit, the magnitude of which will depend on the direction from which Earth is struck. As a result of the nonzero eccentricity, the speed of Earth relative to the Trojan continues to change over the course of an orbit, causing additional instability.

The calculation of the main text shows that direct collisions between incoming rocky bodies and Earth can lead to the disruption of Trojan orbits. In addition to direct impacts, close passages of rocky bodies will also lead to orbital perturbations of Earth. Since these close passages do not add mass to Earth, these encounters will enhance the disruptive effects already considered (i.e., the disruption found in the main text should be considered as a lower limit).

We can estimate the magnitude of this effect. Consider a field of rocky bodies that are streaming toward Earth. As a first approximation, we consider the rocks to have a given mass m, speed v with respect to the Earth, and assume that all of them are coming from the same direction. One could generalize to distributions of (m, v), but we start by comparing the velocity perturbations produced by direct collisions with those induced by close encounters.

The rate at which rocky bodies collide with Earth is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}={nv}\pi {f}_{{\rm{g}}}{R}_{\oplus }^{2},\end{eqnarray} \tag{ C1 }$$

where n is the number density of rocks and where f_g > 1 is the gravitational focusing factor. It will be useful to work in terms of the number of impacts, N, so we write

$$\begin{eqnarray}&&({\rm{\Delta }}t)=\displaystyle \frac{N}{{nv}\pi {f}_{{\rm{g}}}{R}_{\oplus }^{2}}.\end{eqnarray} \tag{ C2 }$$

If all of the collisions have the same sign, the velocity perturbation (to Earth) produced over time (Δt) has the form

$$\begin{eqnarray}&&{({\rm{\Delta }}v)}_{\mathrm{dir}}\approx {{nv}}^{2}\left(\displaystyle \frac{m}{{M}_{\oplus }}\right)\pi {f}_{{\rm{g}}}{R}_{\oplus }^{2}({\rm{\Delta }}t)={Nv}\left(\displaystyle \frac{m}{{M}_{\oplus }}\right).\end{eqnarray} \tag{ C3 }$$

However, we are interested in the square of the velocity perturbation. If the incoming rocks produce both positive and negative angular momentum perturbations, as expected, then the velocity perturbation grows as a random walk (Binney & Tremaine 2008) so that

$$\begin{eqnarray}&&({\rm{\Delta }}v{)}_{\mathrm{dir}}^{2}\approx {{Nv}}^{2}{\left(\displaystyle \frac{m}{{M}_{\oplus }}\right)}^{2}.\end{eqnarray} \tag{ C4 }$$

For completeness, we note that if the rocks all arrived from a single direction, then the expression would pick up an additional factor of N.

Now consider the case of near-misses or close encounters. First we note that while the direct collisions produce velocity perturbations in the direction of motion of the rocky impactor, close encounters produce velocity perturbations in the perpendicular direction. Working in the impulse approximation, the velocity perturbation produced by a fly-by encounter (Binney & Tremaine 2008) has the form

$$\begin{eqnarray}&&\delta v=\displaystyle \frac{2{Gm}}{{bv}},\end{eqnarray} \tag{ C5 }$$

where b is the impact parameter of the interaction (we suppress subscripts, but keep in mind that the direction is perpendicular). For close encounters, the interactions are equally likely to occur on either side of the planet, which means that the velocity perturbations have both positive and negative signs. Moreover, the size of the perturbation depends on the impact parameter. The rate of encounters with impact parameter between b and b + db is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{b}={nv}2\pi {bdb}.\end{eqnarray} \tag{ C6 }$$

As a result, the square of the velocity perturbation accumulates according to the relation

$$\begin{eqnarray}&&d\left[\displaystyle \frac{d}{{dt}}({\rm{\Delta }}v{)}_{\mathrm{cl}}^{2}\right]={\left(\displaystyle \frac{2{Gm}}{{bv}}\right)}^{2}{nv}2\pi {bdb}.\end{eqnarray} \tag{ C7 }$$

After integrating over all impact parameters, and over time, we obtain

$$\begin{eqnarray}&&({\rm{\Delta }}v{)}_{\mathrm{cl}}^{2}=\displaystyle \frac{8\pi {G}^{2}{m}^{2}n({\rm{\Delta }}t)}{v}\mathrm{ln}\,{\rm{\Lambda }},\end{eqnarray} \tag{ C8 }$$

where

$$\begin{eqnarray}&&{\rm{\Lambda }}\equiv \displaystyle \frac{{b}_{\max }}{{b}_{\min }}.\end{eqnarray} \tag{ C9 }$$

The minimum impact parameter must be larger than Earth so that ${b}_{\min }\geqslant {R}_{\oplus }$. In general, we expect b_min $\approx \sqrt{{f}_{{\rm{g}}}}{R}_{\oplus }$. As a working approximation, we can take the maximum impact parameter to be the Hill radius, i.e.,

$$\begin{eqnarray}&&{b}_{\max }\approx {R}_{H}={a}_{\oplus }{\left(\displaystyle \frac{{M}_{\oplus }}{3{M}_{\odot }}\right)}^{1/3}\approx 0.01\,\mathrm{au}.\end{eqnarray} \tag{ C10 }$$

As result, the so-called Coulomb logarithm, $\mathrm{ln}\,{\rm{\Lambda }}\approx 5$. Removing the factor of (Δt) in favor of N, we find

$$\begin{eqnarray}&&({\rm{\Delta }}v{)}_{\mathrm{cl}}^{2}=\displaystyle \frac{8{G}^{2}{m}^{2}N}{{v}^{2}{R}_{\oplus }^{2}{f}_{{\rm{g}}}}\mathrm{ln}\,{\rm{\Lambda }}.\end{eqnarray} \tag{ C11 }$$

Even though the perturbations are in the parallel and perpendicular directions, for direct and close encounter interactions, respectively, we are nonetheless interested in their relative magnitude. Since the direct collisions are also subject to a random walk, the ratio of the two contributions is given by

$$\begin{eqnarray}&&{ \mathcal R }=\displaystyle \frac{({\rm{\Delta }}v{)}_{\mathrm{cl}}^{2}}{({\rm{\Delta }}v{)}_{\mathrm{dir}}^{2}}=\displaystyle \frac{8\mathrm{ln}\,{\rm{\Lambda }}}{{f}_{{\rm{g}}}}{\left(\displaystyle \frac{{{GM}}_{\oplus }}{{R}_{\oplus }{v}^{2}}\right)}^{2}\approx \displaystyle \frac{40}{{f}_{{\rm{g}}}}{\left(\displaystyle \frac{{v}_{\mathrm{esc}}}{v}\right)}^{4}.\end{eqnarray} \tag{ C12 }$$

The distribution of incoming speeds (Figure 1) shows a broad peak near v ∼ 25 km s⁻¹, whereas the escape speed of Earth v_esc ∼ 11 km s⁻¹. The gravitational focusing factor ${f}_{{\rm{g}}}=1+{({v}_{\mathrm{esc}}/v)}^{2}$, so that the ratio ${ \mathcal R }\approx 1$. In other words, in the absence of asymmetries in the collision stream, the velocity perturbations due to close encounters are roughly comparable to those induced by direct collisions. The effects considered in this paper thus represent a lower limit to the full degree of Trojan disruption. Moreover, since the close encounter perturbations are perpendicular to the perturbations produced by direct collisions, the contributions add in quadrature, so we expect the enhancement factor to be $\sim \sqrt{2}$.

This appendix considers another channel for the disruption of Trojan orbits. Given the reservoir of rocky bodies on Earth-crossing orbits that interact with Earth and provide the late veneer, direct collisions with the Trojans can occur. For a given Trojan, the expected number of collisions is given by

$$\begin{eqnarray}&&{N}_{T}={n}_{T}{\sigma }_{T}v({\rm{\Delta }}t),\end{eqnarray} \tag{ D1 }$$

where ${\sigma }_{T}=\pi {R}_{T}^{2}$ is the cross section for the Trojan and n_T is the number density of rocky bodies large enough to disrupt the orbit. Eliminating (Δt) using Equation (C2), the above expression can be written in the form

$$\begin{eqnarray}&&{N}_{T}={N}_{\oplus }\displaystyle \frac{{n}_{T}{\sigma }_{T}}{{n}_{\oplus }{\sigma }_{\oplus }},\end{eqnarray} \tag{ D2 }$$

where n_⊕ is the number density of objects considered in this paper to provide perturbations of Earth's orbit. Note that as the size, R_T , of the Trojan decreases, the cross section, ${\sigma }_{T}\sim {R}_{T}^{2}$, decreases, but the number density of objects, n_T , large enough to disrupt the body increases. Although the size distribution of Earth-crossing bodies is not well known, the two effects roughly compensate. We can thus obtain an approximate result using n_T = n_⊕ and taking the cross section σ_T to be that of an R = 1000 km body. Using these values, along with the size of Earth and the gravitational focusing factor, we find

$$\begin{eqnarray}&&{N}_{T}\approx \displaystyle \frac{{N}_{\oplus }}{50}.\end{eqnarray} \tag{ D3 }$$

Here, N_⊕ is the number of (large) impactors from the simulations of the paper, where we find N_⊕ = 1 − 15, with typical values ∼5–10. This estimate suggests that direct collisions with Trojans are a less likely source of orbital disruption than perturbations to Earth's orbit, by a factor of 5–10.