MEGASIM: Lifetimes and Resonances of Earth Trojan Asteroids -- The Death of Primordial ETAs?

We present an analysis of lifetimes and resonances of Earth Trojan Asteroids (ETAs) in the MEGASIM data set (Yeager&Golovich 2022). Trojan asteroids co-orbit the Sun with a planet but remain bound to the Lagrange points, L4 (60{\deg} leading the planet) or L5 (60{\deg} trailing). In the circular three-body approximation, the stability of a Trojan asteroid depends on the ratio of the host planet mass and the central mass. For the inner planets, the range of stability becomes increasingly small, so perturbations from the planets have made primordial Trojans rare. To date there have been just two ETAs (2010 TK7 and 2020 XL5), several Mars Trojans, and a Venus Trojan discovered. The estimated lifetimes of the known inner system Trojans are less than a million years, suggesting they are interlopers rather than members of a stable and long-lasting population. With the largest ETA n-body simulation to date, we are able to track their survival across a wide initialized parameter space. We find the remaining fraction of ETAs over time is well fit with a stretched exponential function that when extrapolated beyond our simulation run time predicts zero ETAs by 2.33 Gyr. We also show correlations between ETA ejections and the periods of the Milankovitch cycles. Though Earth's orbital dynamics dominate the instabilities of ETAs, we provide evidence that ETA ejections are linked to resonances found in the variation of the orbital elements of many, if not all of the planets.


INTRODUCTION
Asteroids that co-orbit with a planet near its fourth or fifth Lagrange points (hereafter L4 and L5) are called Trojan asteroids 1 . L4 and L5 are stable in the reduced three-body problem for all the planets in the Solar System. Nearly ten thousand are associated with Jupiter's L4 and L5 points as of 2021, and there are expected to be roughly equal 1 km sized Jupiter Trojans to the number in the main asteroid belt (Yoshida & Nakamura 2005). Trojans have been observed co-orbiting with Venus, Earth, Mars, Jupiter, Uranus and Neptune (Connors et al. 2000;Almeida et al. 2009;Gilmore et al. 2010;Connors et al. 2011 To date, there are two known Earth Trojan asteroids (hereafter ETAs): 2010T K 7 and 2020XL 5 (Connors et al. 2000(Connors et al. , 2011Hollis 2022), which both librate far from the Lagrange points within the space between Earth and L3 on trajectories known as 'tadpole orbits'. The two ETAs that have been observed are not in highly stable regimes and are expected to remain ETAs on a time scale of only tens of thousands of years, as determined by numerical simulations (Dvorak et al. 2012;de la Fuente Marcos & de la Fuente Marcos 2021;Santana-Ros et al. 2022). If libration amplitudes remain small ETAs can remain bound to L4 or L5, though if perturbed enough the ETAs can cross near the L3 Lagrange point, and continue to co-orbit but on the opposite side of the Solar system, which is then classified as a horseshoe orbit. Further perturbations will eventually lead to ejection from the co-orbiting regime entirely. In the restricted three body problem, Murray & Dermott (1999) derived analytical constraints for the maximum radial distance an ETA can travel from Earth's orbit and still remain in the tadpole potential well: where r crit is the radial distance of the ETA from the host orbit, and µ 2 = M ⊕ /M . For the Earth, r crit is approximately 0.0011 au. The two known ETAs, 2010T K 7 and 2020XL 5 , have heliocentric semi-major axis determined to be a = 0.999273055 ± 8 × 10 −9 AU and 1.0009 ± 0.0002 AU, respectively (de la Fuente Marcos & de la Fuente Marcos 2021). Both of these semimajor axis currently are within 1 ± r crit criteria for the Earth, but even these are stable for a short period of time compared with the age of the Solar System.
There are a few possible sourcing mechanisms for ETAs. They could have formed in the proto-planetary disk and fallen into the eventual L4 or L5 (Montesinos et al. 2020). At a slightly later stage, they may have been captured debris produced from major planetary impacts or planetary migration. Finally, they may be sourced from perturbed populations of asteroids elsewhere in the Solar System such as ejected asteroids from the main asteroid belt. The first two options would have occurred early in the formation of the Solar System, and we consider these to be "primordial" in nature. Primordial sourcing through the impact hypothesis may also explain lunar formation. The giant impact hypothesis, where the Earth was impacted by another early planet that came from a 1 AU orbit in the solar nebula has been proposed. Belbruno & Gott (2005) suggest a Mars-sized proto-planet could have formed in a stable orbit among debris at either Earth's L4 or L5 Lagrange point and collided with Earth to form the eventual Earth-Moon system. It is reasonable to expect a cloud of debris from such a collision may have provided a primordial source of ETAs (Canup 2012). If found, such a population of asteroids would offer a unique lens into the evolution of the early Solar System and the Earth and Moon system. This fact has driven much interest in their discovery.
Modeling of the ETA population with n-body simulations first appeared in the literature by Szebehely (1967), with results showing there to be regions where long term stability of ETAs can exist. Since then a plethora of numerical simulations have be carried out, but our focus is on a recent simulations that study the lifetimes of ETAs. JeongAhn & Malhotra (2010) modeled the lifetimes and stability of three distinct populations of ETAs, built by allowing only one of the semimajor axis, eccentricity, or inclination to vary from that of Earth's orbit. These populations were defined as the following: 'A', whose semi-major axes were varied between 0.988 and 1.012 AU, 'E' with eccentricities ranging between 0 and 0.4, and the 'I' population with inclinations varied between 0 and 180°. The 'A' population of ETAs was found to retain 80% of the initial population by 100 Myr. The 'E' and 'I' ETA populations were found to be much more unstable as a whole, with the 'E' population only retaining 20% of the ETAs by 100 Myr, and the 'I' population falling to less than 10%. The 'I' population, unlike the other two populations, was also found to lose approximately half of the ETAs suddenly after 1 Myr of simulation time.
A potential population of ETAs is suggested byĆuk et al. (2012), who analyzed the long term stability of both tadpole and horseshoe orbits in a sample of 1000 asteroid orbits. The initial semi-major axis are allowed to be one of 25 values evenly spaced from 1 to 1.01 AU with 40 possible inclinations evenly spaced from 0 to 40 degrees. They found a fraction of their initial set of asteroids remained stable as Earth horseshoes and tadpole orbits at the end of their 700 Myr Earth simulation. The orbits that remained stable until 1 Gyr were found to exist with semi-major axis between 1 and 1.002 AU, eccentricities between 0 and 0.2 and inclinations from 0°t o 19°.
Marzari & Scholl (2013) applied a frequency map analysis (e.g., Marzari et al. 2002) to an ETA simulation which indicated some ETA orbits are Gyr-stable . In this case, the prediction of long term stability is gleaned from a 10 5 year integration. A frequency map analysis in Keplerian orbital element parameter space was applied to a 12 Myr data set (Zhou et al. 2019). Two stability regions were found for inclinations below 15°and another 24°to 37°. The orbital instabilities for inclinations between 15°and 24°are believed to arise from the ν 3 and ν 4 secular resonances. Nodal secular resonance ν 13 is able to change initial inclinations by up to 10°and the ν 14 resonance could increase inclination 20°or more. Zhou et al. (2019) found that tadpole orbits experience great differences in the frequency spaces of their secular precession and that secular precession in tadpole orbits are more sensitive to the frequency drift of the inner planets. This could indicate a clearance of all tadpole orbits. Furthermore, an asymmetry was found between the prograde and retrograde rotating ETAs possibly due to the Yarkovsky effect (for a historical reference see Beekman 2005) -a thermal effect in which the light from the Sun heats up the asteroids, the re-radiated light causes a small perturbation radially, which dependent upon spin, will result in an inward or outward migration relative to the Sun. This effect produces an increase in orbital instability for small asteroids. If an ETA is to survive 1 Gyr they find it must have a radius larger than 90 m for a prograde spin and 130 m for a retrograde spin. For an ETA to survive the age of the Solar System a radius of 278 m for a prograde spin and 400 m for a retrograde spin is required.
Recently, Christou & Georgakarakos (2021) found similar results toĆuk et al. (2012); Zhou et al. (2019) insofar as the stability of horseshoes is found to be greater than tadpoles, and that horseshoes and tadpoles may easily survive the age of the Solar System. From the simulation data a constraint of 3-12 ETAs with no more than 60 may exist bound to the L4 and L5 points. The simulated tad poles were found to have excitation in the orbital eccentricities that persist until the eccentricity is beyond 0.1 or the semi-major axis deviates from that of Earth by more than 0.0028 au. Horseshoes were instead found to escape through an excitation of the orbital inclination for angles greater than 13°. The simulations provide further evidence to that in Dvorak et al. (2012); Zhou et al. (2019), that secular resonances are responsible for destabilising ETAs. The fraction of ETAs that remain over the simulation appears to fall roughly linearly from 0 to 750 Myr, flattens until 1.3 Gyr before continuing down at a slightly shallower linear rate until the end of the 2 Gyr simulation.
Large N-body simulations that map the full picture of ETA stability and trajectories are important in placing limits on observational attempts to detect ETAs. Wiegert et al. (2000) found the densest location of the ETA distribution to be slightly inside the orbit of Earth near L4 and determined observability of their simulated ETA populations. Todd et al. (2012) further explored optimal search strategies. A case was made that 70-80 ETAs are missing from observations (Christou 2019).
Observers have searched for such primordial populations suggested by the simulations. Lifset et al. (2021) used DECam data and a novel shift-and-stack method to set a tight upper limit on ETAs with H<21.77. They find their observations support an ETA population at L4 of N ETA < 1 for H = 14.28, N ETA < 7 for H = 16, and N ETA < 642 for H = 22. They also reviewed and updated the upper limits of all other known searches for this population, at L4 (Yoshikawa et al. 2018) and the similar L5 (Whiteley & Tholen 1998;Cambioni et al. 2018;Markwardt et al. 2020). Though these results eliminate the chances of a large population of 1 kmclass ETAs at L4, they do not eliminate a population of ∼100 m ETAs, which remain a target of deep interest for Solar System formation and planetary defense.
Here we present results of the highest fidelity simulations of ETAs to date (Yeager & Golovich 2022). The MEGASIM offers a complete coverage of the initial Ke-plerian orbital element parameter space, where as previous work generally fixed multiple parameters to be equal to Earth, the MEGASIM requires no such conditions on the 6-dimensions required to initialize ETA orbits. The analysis in this paper centers on the lifetimes (survivability) with a specific interest in answering whether there can exist a primordial population of such objects. In §2 we describe our model. In §3 we present our analysis and results, and in §4 we offer a discussion of our results compared with the literature and our conclusions.

THE MODEL
We utilize two different integrators to build two distinct data sets. One set of simulations was carried out using the IAS15 n-body integrator and another using the Wisdom-Holman symplectic integrator referred to from here on as WHfast (Wisdom & Holman 1991;Rein & Tamayo 2015). Both integrators are a part of the REBOUND package (Rein & Liu 2012;Rein & Spiegel 2015). The Solar Systems are initialized with the Sun and eight planets, Mercury, Venus, Earth-Moon, Mars, Jupiter, Saturn, Uranus and Neptune. The initial planetary positions are gathered from the JPL HORIZONS system at an epoch of JD2458475.30035 (Giorgini et al. 1996). We utilized the Catalyst, Quartz and Ruby clusters managed by Livermore Computing at Lawrence Livermore National Laboratory 2 .
The two ETA data sets in this paper will be referred to by the integrator that was used (IAS15 and WHfast, respectively). Initially only the IAS15 integrator was used to generate the orbits of ETAs and was chosen because of its ability to use a variable time step meaning near approaches to Earth would be well-resolved. While the IAS15 integrator lends to very high accuracy it comes with a high computational cost, limiting the duration simulations could reach. The choice was then made to carry out another set of simulations with the WHfast integrator to reach much longer time scales. The REBOUND WHfast integrator takes several parameters to define its use, a time step of 1 day, symplectic correctors were set to 17th order (sim.ri whfast.corrector = 17) and the sim.ri whfast.safe mode parameter was turned off. Runs were initially conducted without the use of sympletic correctors, and it was found that the stable orbits of the ETAs was poorly resolved. This comes with a benefit of being able to compare the results from both integrators using large number statistics to better understand any possible limitations of using the WHFast integrator.
Parallelization of the simulation is achieved by running 44,800 initializing batches of 500 asteroids injected with six orbital parameters selected from a random distribution, as described below. The number 500 was chosen to balance running enough asteroids so that individual simulations would have asteroids in them without wasting CPU time integrating mostly planets and running few enough asteroids so that the 24 hr wall clock limit for the high-performance computing (HPC) clusters was not prohibitive of simulation progress. Each model is then integrated forward with REBOUND. Position and velocity information for each is output every 1000 years but the integration time steps are shorter. The IAS15 simulations were carried out until either no stable ETAs remained or until a simulation of 50 million years was achieved. The set of simulations integrated with WHFast are not cut off until they either have lost all stable ETAs or until a simulation time of 1 Gyr was reached.
Removal of ETAs occurs if they drift beyond the Earth -L3 line, effectively removing any horseshoe or nonco-orbiting objects. Asteroids that were ejected due to complete instability eventually crossed this line and were removed.
The ETA orbits are initialized with a true longitude (θ) randomly and uniformly distributed between Earth and L3 on the side of Earth's orbit containing L4. The argument of pericenter and longitude of ascending node are randomly assigned values between 0 and 2π assuming a uniform distribution. The semi-major axis (a), eccentricity (e > 0) and inclination (i) were sampled from a three dimensional Gaussian with mean and covariance (2) A visualization of the 11.2 million initialized orbits is provided in Figure 1. The initial positions of the inner planets is also provided, with circular pseudo-orbits overlaid for reference.
The initial Keplerian orbital elements extended well beyond ranges expected for ETAs. Importantly, none of the orbital elements are fixed to that of Earth, which physically may not be required for stability. This is a markedly different setup from the simulations in the literature, which all fixed a subset of the orbital elements and often varied only one a time. While our initial parameters included many orbits we are not interested in, our method for removing asteroids from the simulation does not allow any Earth co-orbits that cross the Earth-Sun-L3 line to remain in the simulation for very long. For a horseshoe orbit to exist within the simulation, it would have to find itself on the L4 side of the Earth-Sun-L3 line every 1000 years to avoid removal from the simulation. This is exceedingly unlikely, and we took care to remove those that did survive a few timesteps by chance from our results. Our broad initialization does result in a computationally inefficient simulation, but in the end it results in the broadest look at ETA stability at high statistical fidelity over a Gyr, which is truly unprecedented in the literature. Figure 2 shows a histogram of binned ejections over time for the first 50 Myr of each simulation. The loss of asteroids is very similar for both integrators. The height of the bins corresponds to the number of ETAs that are removed from the simulation with bin widths of one hundred thousand years. There is a periodic oscillation in the number of ETAs removed with time, apparent from the height of the bins with a period of roughly 2 Myr. This finding is discussed further in §3.2.

Comparing the two simulations
The simulations begin with 11.2 million ETAs and after every simulation time step (1000 yr) the number of ETAs remaining in the simulation is recorded. The resulting curve, shown in panel (a) of Figure 3, is the number of remaining ETAs counted every 1000 years divided by the number initialized. An ETA lifetime is then defined as the amount of time the ETA spent in the simulation before being removed due to crossing the plane perpendicular to the ecliptic that contains the Earth and L3; that is, asteroids co-orbiting with Earth but leading by no more than 180°.
In the early stages of the simulation the fraction of remaining ETAs rapidly decreases to fewer than 10% within the first ∼10 kyr. This initial phase of ejections is indicative of our broad initialization. We earlier drew a distinction between these loosely bound interlopers 3 and the long term ETAs that we seek to study here. After 10 kyr, the ejection rate slows for the majority of the simulation.
The absolute difference between the two integrators (IAS15 and WHfast) remaining ETA fractions is provided in panel (b) of Figure 3. The typical difference between the two integrated data sets ETA fractions is of order 10 −5 , which corresponds to ∼100 out of 11.2 million ETA particle difference between the sims. We initially worried that the use of the fixed time step integrator for WHFast would lead to significant differences in 3 There is value in studying these objects as they closely resemble the two known ETAs. Future work will study this early ejected population.  long term stability estimation as the effect of near Earth approaches may have been smoothed over. These results show that the WHfast integrator works extremely well at matching the results of the much more computationally intensive IAS15 integrator for timescales of 50 Myr. The deviation between the two models remains very consistent even at 50 Myr indicating the accuracy is maintained on significantly longer time scales. The presence of the periodic ejections in both panels of Figure 2 further suggests that the two simulations are capturing the same physics to high fidelity. Upon finding this, we cut off integration for the IAS15 simulation at 50 Myr and carried forward the WHFast simulation alone. The remainder of our results will focus on the WHFast results.

Short Term Ejection Signatures
The periodic ejection pattern apparent in 2 has a large amplitude. Ejections vary by as much as a factor of ten with a rough period of ∼2 Myr. Panel (a) of Figure 4 shows this continues for WHFast to beyond 100 Myr. At 150 Myr the periodicity begins to wash out and by 200 Myr it no longer apparent. Panel (b) of Figure 4 shows this periodicity at higher time-resolution with 1000 year bins. There are higher frequency oscillations apparent as well indicating a blend of ejection sources.
Here we explore the evolution of the individual parallel simulations relative to one another. In order to simulate so many asteroids, we initialized 22,400 identical Solar Systems each with 500 ETAs. Figure 5 shows how the position of the planets in each simulation deviate from one another over time. The 100-200 Myr time-frame matches well with the diminishing of the periodic ejection pattern apparent in Figure 2. This demonstrates that the oscillatory behavior is likely always present; however, it is the coherence of the individual Solar Systems diminishing rather than the periodicity to cause the cumulative binned ejections to change appearance beyond 150 Myr.
In Figure 6 we present a Fast Fourier Transform (FFT) of the ETA ejections with a number of important frequencies. The FFT is computed from the bin counts of a histogram of ETA lifetimes with bin widths of 1000 years. In the top panel of Figure 6, we show the alignment of several peaks in the 65 kyr to 1 Myr range with the well-known Milankovitch cycles, which were first considered to describe changes in Earth's climate and associated with variations in the orbit of Earth. The periods for known changes to Earth's eccentricity, inclination and pericenter are overlaid on the ETA ejection FFT. Periods due to changes in eccentricity are shown as orange vertical lines at 95, 125 and 413 kyr, which loosely combine to a period of 100 kyr, green lines are periods of inclination changes at 70 and 100 kyr, and the purple line is the period of apsidal precession at 112 kyr (Muller & MacDonald 1997;Laskar et al. 2011). It is not surprising that Earth's orbital dynamics are a driving factor in ETA ejections, but to our knowledge this is the first indication that the same variations in Earth's orbit can be linked to both climate dynamics and ETA dynamics.
In the bottom panel of Figure 6, we compare the ejection FFT with s and g resonances of each planet. The periods of the planetary s (solid) and g (dashed) resonances overlaid on the ETA ejection FFT, spanning 10 kyr to 10 Myr are colored by planet. The resonances for the inner planets are from Laskar (1990) and the outer planets Nobili et al. (1989). The s and g resonances of the planets apart from Uranus do not directly match up with any peaks in the ETA ejection FFT. The resonances of Neptune may play a role in the ETA ejections as they are in the range of the multiple ETA ejection peaks centered around 2.3 Myr that we discussed above. The s resonance of Jupiter is found at 129 Myr off the right of panel (b), but there is no correlated peak in the ETA ejection FFT, so we truncated the figure.
Next we compare the FFTs of all orbital elements for each planet and the ETA ejections. The planet orbital element FFTs are computed using a single Solar System simulation spanning 1 Myr to 1 Gyr. We checked for consistency across a handful of our parallel Solar Systems and found good qualitative agreement. Figure 7 shows the FFT for six Keplerian orbital elements for each planet (separated by subplot) along with the FFT of the ETA ejections for comparison. The vertical axis labels indicate the Keplerian orbital element, but for further guidance, they are from top to bottom: longitude of the ascending node (brown), argument of periapsis (purple), true longitude (red), inclination (green), eccentricity (orange) and semi-major axis (blue). At the bottom of each panel, we show the ETA ejection FFT for comparison (black). The vertical axis is scaled such the the highest power for a given FFT is of equal height -i.e., the relative heights between the various FFTs should not be compared.
The Milankovitch cycles shown in Figure 6 appear in the orbital FFTs of Earth's orbital elements. Earth's eccentricity shows the three modes near periods of 95, 125 and 413 kyr and the inclination at 70 kyr, though the 100 kyr period is a very small peak. The simulation orbital element data FFTs correspond well to the spread of frequencies seen in ejection FFT (black). However, the peak with a period just larger than 400 kyr lines up not only with Earth's eccentricity mode, but also strong signal in Venus's eccentricity FFT as well as Uranus's and Neptune's inclination FFT. Figure 8 is presented analogously to Figure 7 for FFT signals with longer periods. The longest period of ejection oscillation peaks near 2.3 Myr with the gradual increase in FFT power from 1.5 Myr to 2.75 Myr. This ∼ 2 Myr period corresponds to the oscillation seen in Figure 2 and 4. The full structure of the ETA ejection FFT seems to be built from an additive combination of many planetary resonances, involving all components of planetary orbital variation.

Long Term Ejection Signatures
Previous work found fall off rates of initialized ETAs of various functional forms including logarithmic and and   Our simulation is different in that we initialized a much broader region of orbital space, and we simulated many more asteroids out to 1 Gyr. In the literature, the ETAs were initialized on elliptical orbits with precise Earth-like properties. For example, in previous simulations the longitude of the ascending node and argument of perihelion have not been varied from that of Earth. Initializing ETA orbits over the full Keplerian parameter space provides a full accounting of the zones of stability over long time-scales, the orbital elements from which they are possibly initialized, and their associated lifetimes. It is our primary goal to study this broader population of potential ETAs.
We consider here the fractional remaining population as a function of time at the time resolution of our simulation output (1000 year increments). For the first ∼10 time steps, we find a logarithmic fall off as the most unstable initialized orbits are quickly cleared out (see the leftmost portion of the black curve in the top panel of Figure 3). By 15 kyr the rate of change in the surviving fraction of ETAs departs from the rapid logarithmic fall off and is well modeled by a stretched exponential from 15 kyr to 1 Gyr. This is a function of the form where A, λ, β and C are free parameters. We show the fit with a red dashed curve in Figure 9. The black line is the WHFast simulation data with the thick black portion indicating the time span for which the function was fit. The fit, when extrapolated beyond the simulation time predicts zero asteroids (we will call this time t 0 ) by t 0 = 2.33 Gyr. We will discuss the implications of this in §4.
To assess the uncertainty in our stretched exponential fit and extrapolated time to zero remaining asteroids, we preformed a bootstrap re-sampling analysis on the asteroid lifetimes. For each initialized asteroid, we stored the time it was either ejected or 1 Gyr for those that survive the duration. To create bootstrap samples of the data, we selected 11.2 million asteroid lifetimes from our data (with replacement) and fit a stretched exponential to the realized fractional remainder as a function of time. For the full data set, our bootstrap analysis gives a mean and 80% confidence interval of t 0 = 2.33 +3.04×10 −3 −4.55×10 −4 Gyr. Thus our simulation provides evidence that a primor-dial population of ETAs that were leftover from the Solar System formation (i.e., that were around when the planets assumed their current arrangement as our simulation was initialized) would not survive to the current age of the Solar System.
Since our extrapolated lifetime is less than the age of the Solar System, we performed an analysis on the estimated zero time (t 0 ) value of a stretched exponential fit as a function of the amount of our simulation data used. That is, we fit a stretched exponential to our bootstrap samples that are truncated at 100 Myr intervals out to 1 Gyr. For all fits, we have fit from 15 kyr to the truncated time. The results of this procedure are presented in Figure 10. The black dot indicates the mean t 0 (vertical axis) predicted for all curves fit to a given range of data, and the orange bars span the 80% confidence interval for the bootstrap samples. As more data is fit, the stretched exponential tightens to a t 0 under the age of the Solar System. We will discuss this further in §4. The mean and 80% confidence intervals for all parameters of Equation 3 and t 0 for the bootstrap analysis are presented in Table 1. Figure 7. FFT of planetary orbital elements using 1 Myr to 1 Gyr of simulation data. Panels only span 10 kyr to 1 Myr to focus on the many lower frequencies that exist in both the planetary orbital elements and ETA ejections. In each panel is a different planet with line colors from bottom up: black/ETA ejections, blue/semi-major axis, orange/eccentricity, green/inclination, red/true longitude, purple/argument of periapsis, brown/longitude of the ascending node.  (3) can be physically interpreted as one over the averaged relaxation time of the relaxing species, in our case the entire ETA population. Based on the best fit β values of Table 1, a corresponding relaxation time of the ETAs increases with time, ranging from 6.6 Gyr with a 100 Myr data fit up to 8.8 Gyr when fitting to the full 1 Gyr of data. Such a long relaxation time in the ETA population would indicate that chaos will eventually lead to a depletion of all ETAs, as our stretched-exponential fit suggests. The perturbations to the ETA orbits cause both short and long timescale diffusion of the ETAs. Thus they can never truly reach an equilibrium before further perturbations occur. Figure 9. The remaining ETA fraction as a function of simulation time in black with a best fit stretched exponential (red dashed line) overlaid. y(t) and t refer to variables found in Equation (3). The area where the black line is thicker is used for the stretched exponential fit to avoid the initial stages of the simulation where very unstable ETAs filter out. The fit spans 15 kyr to 1 Gyr.

The physics of the stretched exponential
The number of ETAs remaining with time is found to follow a stretched exponential form. Here we explore the literature of stretched exponential functions, which have been used frequently in scenarios similar to ours. Dobrovolskis et al. (2007) fit the survival of simulated particles ejected from the Saturnian Moon Hyperion with a stretched exponential, and also compiled a list of many other related astrophysics particle removal simulations where the stretched exponential appears (Gladman & Duncan 1990;Holman & Wisdom 1993;Levison & Duncan 1994;Gladman et al. 1995;Dones et al. 1996;Gladman et al. 1996;Farinella et al. 1997;Gladman 1997;Holman 1997;Burns & Gladman 1998;Dones et al. 1999; Figure 10. The range of predicted by the stretched exponential fits to bootstrap resampled data as the domain over which the fit is performed is increased. The fits all begin at 15 kyr and extend to 100 to 1 Gyr incremented by 100 Myr. The width of the orange lines extends from the 10th to 90th percentile obtained from bootstrap samples of remaining ETA distribution. Evans & Tabachnik 1999;Malyshkin & Tremaine 1999;Dobrovolskis et al. 2000;Gladman et al. 2000;Hartmann et al. 2000;Tabachnik & Evans 2000;Robutel & Laskar 2001;Armstrong et al. 2002;Chambers et al. 2002;Nesvornỳ & Dones 2002;Quintana et al. 2002;Nesvornỳ et al. 2003;Dobrovolskis & Lissauer 2004;Alvarellos et al. 2005;Gladman et al. 2005;Zeehandelaar & Hamilton 2005;Hamilton & Zeehandelaar 2006;Zahnle et al. 2007;Alvarellos et al. 2008).
Our further literature review found stretched exponential used to model physical processes in many different areas of research, with the common thread being underlying chaotic dynamics. We provide here an incomplete but expansive list to demonstrate the breadth of the literature: • General disordered systems (Love 2019) • Li 7 nuclear spin-lattice relaxation (Trinkl et al. 2000;Kaps et al. 2001;Johnston et al. 2005) • Chaos in the formation of cosmological large scale structure  • Cosmological decay of Higgs-like scalars (Boyanovsky & Herring 2019) • Thermal relaxation of self gravitating sheet model (Joyce & Worrakitpoonpon 2010) • Solar flare distributions (Najafi et al. 2020) • Star-forming cloud models (Collins et al. 2012) • Sand drift wind models (Yizhaq et al. 2020) Table 1. The mean and 80 % confidence intervals for the stretched exponential parameters and the extrapolated time to zero remaining ETAs. t end gives the truncated time over which the fits were performed. All fits begin at 15 kyr. A, λ, β, and C are the parameters of the stretched exponential of form: Equation (3). t0 is when each curve crosses 0 on the y-axis. Values provided are the mean of the bootstrap population with their 80% confidence intervals. The confidence intervals span the 10-90 quantile range for the bootstrap samples of our ETA simulated lifetime data.
• Viscous remnant magnetization dating method applied to tsunami boulders (Sato et al. 2019) • Buoyancy driven turbulence (Bershadskii 2016) • Distributed chaos and isotropic turbulence (Bershadskii 2015) • Modeling of earthquake aftershocks (Gasperini & Lolli 2009;Mignan 2015) • Rainfall distribution (Wilson & Toumi 2005) • Solid-state physics for creep, annealing, electrical impedance, and magnetic spin relaxation (Palmer et al. 1984;Peterson 1989;Halsey & Leibig 1991;Scher 1991;Phillips 1996) Based on the literature review above, it is not surprising that our ETA remainder fraction is well fit by a stretched exponential, but what is the origin of the mathematical form?
Here we draw a connection between our findings and a simplified model of the diffusion processes. Dobrovolskis et al. (2007) suggested diffusion in the space of orbital elements was related to their findings regarding the ejected material from Hyperion. Simplifying from the stretched exponential, we first consider the onedimensional diffusion equation: where P (x, t) is distribution of the diffusing fluid or particles. D is the diffusion constant and incorporates the physics involved with a given diffusion process. This can be separated into spatial and time components (X and T , respectively): which have solutions: respectively, where a, b, and c are constants related to the boundary conditions, and k is the separation constant. The time-dependence is exponential in nature, with a characteristic relaxation time related to the diffusion constant D and the separation constant k.
Continuing with the idea that the ejection of ETAs is analogous to diffusion processes, let us assume some function f (t) is equal to a series of exponential functions: where C n and a n are for now just some constants for fitting. The philosophy here is that, following Dobrovolskis et al. (2007), the diffusion processes in the space of orbital elements was related to the stretched exponential functional form. We know that orbital elements in the many-body gravitational problem can induce resonances that perturb minor planets (e.g., the Kirkwood Gaps in the main asteroid belt Kirkwood 1891). However, there are also subtler resonances that occur at non-regular intervals that perturb the ETA population with individual diffusion processes. The derivatives of this formulation is then: C n a n n exp{a n t} (10) The matrix form of this equation is: where the square matrix is known as the Vandermonde matrix, which itself has many interesting connections in applied mathematics. Relevant for our discussion here is the connection between the Vandermonde matrix and the discrete Fourier transform matrix (Bäckström et al. 2014). Furthermore, the Vandermonde matrix appears in the theory of Dyson Brownian motion (Eichelsbacher & König 2008) 4 . The relationship to Brownian motion here is especially interesting as it implies the diffusion of ETAs out of the loosely bound potential is associated with a random walk due to accumulating perturbations. These perturbations are caused by consecutive resonances with the planets over long timescales. Brownian motion is often invoked as a random walk in a flat potential, but a potential is not excluded mathematically, and so we interpret our ETA ejection process as each ETA experiencing an accumulation of random perturbations, which are linked with independent diffusion processes with their own characteristic relaxation time. The ETAs in our initial population are broad, and often easily ejected from a random accumulation of these perturbations. Others withstand the random perturbations for long timescales; however, over time the accumulation tends to all the ETAs by a finite timescale. It has also been shown that a stretched exponential can be represented as a continuous spectrum of pure exponential functions (Johnston 2006): where s ≡ λ * /λ and λ * is the characteristic relaxation time for a given pure exponential and P (s, β) is a probability distribution such that ∞ 0 P (s, β) ds = 1. In this formulation, the parameter β of Equation (3) is related to the logarithmic full-width at half maximum (FWHM) of the probability distribution P (s, β). P (s, β) has a cutoff at low s that decreases monotonically with decreasing β such that β can measure the small relaxation-rate cutoff of P (s, β).
Our finding that there are no simple discrete set of resonances from the planets in Figure 6, 7, and 8 to describe the complicated nature of the ETA ejection FFT further bolsters our discussion here, where we tie the ETA ejections to Brownian motion. The random walk nature of Brownian motion is analogous in some sense here to the accumulation of random perturbations experienced by the ETAs over time. As we showed in §3.2.2, the β value from our stretched-exponential fit suggests they can never reach equilibrium from prior perturbations before future ones occur.

The Yarkovsky Effect
Because we do not account for additional forces beyond n-body gravity the lifetimes of the ETA orbits may be considered as upper bounds for stability. This further strengthens the argument for no existing primordial population of ETAs today. The strength of the Yarkovsky effect and thermal radiation from the Sun depends on the rotation speed, size, albedo and density of a given asteroid. Zhou et al. (2019) found smaller asteroids with radii below 90 m or 130 m (depending on if the orbit is rotating prograde or retrograde, respectively) are cleared out of the ETA regime on timescales less then 1 Gyr. This finding suggests that any primordial population of ETAs is likely composed of asteroids larger than 100 meters.

The primordial population of ETAs
We would like to point out here that we are particularly focused on primordial ETAs, but it is worth exploring exactly what that means. While it is easy to invoke ETAs coming from the early Solar System, some consideration is needed to think through the implications for the results in this paper.
Since we initialized our asteroids with Earth-like orbits in a Solar System with present day planets, we either have implicitly assumed that asteroids survived migration of the gas giants or that an additional mechanism for sourcing these asteroids occurs after migration. It certainly seems reasonable to assume that very early Solar System ETAs may have been disrupted by gas-giant migration since our findings show that much smaller perturbations eventually can clear out ETAs anyway, so it is a reasonable step to assume much more rapid large scale motion by the planets likely had an out-sized effect on the ETA population.
A single model of the Solar System that can account for each of gas-giant formation, Jupiter Trojan asteroids, the asteroid belt, the Kuiper Belt, the present arrangement of the inner and outer planets, as well as the observed impact histories of the inner planets and the Moon does not yet exist. Several models have succeeded in modeling a subset of the above. One leading model, colloquially named the Nice model (after the Southern French city) provides an explanation for the arrangement of the gas-giants, a means of stripping the protoplanetary disk of its gas in order to generate rocky planets , the distribution of asteroids from the asteroid belt to the Kuiper belt , and an explanation for the apparent late heavy bombardment (LHB) 3.8 − 4.1 Gyr ago from lunar samples . There are other models, and even extensions to the Nice model, but we focus on this one because of last point above. The fact that evidence exists of massive impacts of large asteroids ∼4 Gyr ago in the inner Solar System after the formation of the Earth-Moon system suggests another potential sourcing term for ETAs. Gomes et al. (2005) suggests impacts to the moon by asteroids as massive as 8 × 10 21 g. Assuming 2 g cm −3 , which is suggested in Figure  7 of Carry (2012) for asteroids of this mass, this is an asteroid of ∼ 200 km in diameter. The escape velocity of the moon is ∼ 2.4 km s −1 , so it is feasible that LHB sourced inner Solar System Trojan asteroids including ETAs through impacts such as those suggested here, but also with other inner planets, moons, or minor planets.
We still stand by our assessment that if ETAs were sourced by LHB or earlier from the giant impact hypothesis or any other early Solar System mechanisms that were able to survive to the present arrangement of the Solar System, these types of ETAs would be considered "primordial" either way. Each of these occurred before the Solar System was 1 Gyr old leaving as much as 3.5 Gyr to survive to the present day. While previous work has shown ETA stability on the age of the Solar System is possible (Ćuk et al. 2012;Zhou et al. 2019), our simulations point toward it being unlikely.
Of the 11.2 million orbits initialized in the WHfast data set, 22,441 ETAs total survive a 1 Gyr integration. We showed in Figure 10 that increased integration extrapolated to more precise lifetimes shorter than the age of the Solar System (see Table 1). Extrapolation of the stretched exponential obtained from the full 1 Gyr of data predicts zero ETAs will remain in the simulation by 2.33 Gyr.
To compare, there are roughly 20,000 near-Earth objects (NEOs) larger than 100m in the inner Solar Sys-tem (Mainzer et al. 2011(Mainzer et al. , 2014. This makes a good comparison for our surviving ETAs, since in §4.2 we argued that our simulation should be interpreted only to contain those larger than 100 meters because anything smaller would have been even more rapidly ejected due to the Yarkovsky effect. Primordial ETAs are unlikely to be more numerous than NEOs of the same size, which have a constant source in asteroids ejected from the main asteroid belt. Thus, we argue that zero ETAs in our simulation is indicative of zero primordial ETAs today 5 .

Extrapolating the stretched exponential fit
As we continue to run the simulations, we will be able to put the predicted zero time to the test by comparing the actual time when the simulations run out of asteroids to the predicted zero time. We note that should all of the simulations run out of asteroids, it does not directly indicate that the Solar System would run out of asteroids, but rather that with the current set of simulations the remaining asteroid fraction today would be less than one in 11.2 million of the initial population. This raises the question: is 11.2 million enough? Indeed, 11.2 million was picked for our own practical reasons more than anything. It is a large number, larger than the literature has seen by a wide margin, but also it fit on a manageable amount of nodes on our HPC to run them in entire batches to completion within our HPC wall-clock limit. In principle, one could simulation the settling of debris into the L4 or L5 points in a proto-planetary disk, or the settling of debris from later stage large impacts, but this is beyond the scope of our work, and we know of no such simulation in the literature.
Finally, the zero time for the best fit stretched exponential (see Equation 3) crosses zero because the constant term C is negative. To obtain a fit at timescales of 10 kyr and 1 Gyr the constant factor is required. Without it equation 3 does not fit the simulation data well at both 10 kyr and 1 Gyr timescales, nor does it capture the steep downturn in remaining ETAs leading up to 1 Gyr. Without an additive constant an exponential has an asymptote on the x-axis never reaching zero. Which if taken at face value would imply an exponentially decreasing population would survive over arbitrarily long timescales. In practice a cut off should be made on the remaining ETA fraction, as below some fraction the likelihood of at least 1 remaining ETA becomes effectively 5 We will explore the potential of ground based surveys for both the surviving primordial ETA population as well as interloping ETAs like 2010T K 7 and 2020XL 5 , for which we have collected a lot of information on in MEGASIM -namely those that eject in time scales similar to the estimated lifetime of the known interloping ETAs.
zero. Where such a cut off should be made requires a constraint on the primordial mass that was initially bound to L4. As we mentioned in §4.3 the primordial L4 mass is not well constrained by Solar System models.

Conclusions
HPC at Lawrence Livermore National Laboratory has provided us the opportunity to build the largest ever simulation of ETAs initialized over the full Keplerian orbital parameter space with unprecedented fidelity. In this work we presented an analysis of the population of ETAs in regards to their lifetime as a population. Using the WHfast integrator the orbital data sets span an integration time of 1 Gyr for an initialized population of 11.2 million asteroids and all eight planets. We have continued integrating the simulation, but our results here stand for themselves. Further analysis is planned for our simulated data set. Coming papers will cover observability in current and future astronomical surveys such as the Zwicky Transient Facility and the Vera C. Rubin Observatory's Legacy Survey of Space and Time. We will also explore Earth impact likelihoods and the trajectories of quasi-stable interlopers like the two known ETAs.
Here we conclude our discussion with a brief list of the main findings of our analysis in the context of our literature review.
• We predict no ETAs can survive to the present day if those asteroids were initialized in the early Solar System.
• The fractional remainder of ETAs in our simulation over time follows a stretched exponential form. Extrapolating the remaining population of ETAs predicts zero by 2.33 Gyr.
• We performed an extensive review of stretch exponential functions in the literature, and we find similarities in those cases to ours here, and we related the physics of stretched exponentials to additive diffusion processes that are ultimately unable to relax in time before they are all ejected.
• Using a FFT analysis of the ejected ETAs and the planet's orbital dynamics, we showed that portions of the ejection pattern of ETAs aligns well with changes in Earth's orbital elements that cause the Milankovitch cycles that have been indicated in the climate record of Earth. The most important orbital variations seem to be in the semi-major axis, eccentricity and inclination of Earth.
• We showed that broad signals in the ejection pattern of the ETAs align with many resonances in all of the planet's orbital elements, which indicates the complexity of the physical process and also leads to the stretched exponential functional form and our interpretation.
• We continue to integrate the simulations beyond 1 Gyr. When all asteroids have been ejected or if the functional form of the surviving ETAs changes our interpretation, we will follow up on our analysis.