# MEGASIM: Distribution and Detection of Earth Trojan Asteroids

Travis Yeager Lawrence Livermore National Laboratory
7000 East Ave, Livermore
Livermore, CA 94550, USA  Nathan Golovich Lawrence Livermore National
Laboratory
7000 East Ave, Livermore
Livermore, CA 94550, USA

###### Abstract

Using the MEGASIM, we present spatial distributions of Earth Trojan Asteroids
and assess the detectability of the population in current and next-generation
ground-based astronomical surveys. Our high-fidelity Earth Trojan Asteroid
(ETA) distribution maps show never-before-seen high-resolution spatial
features that evolve over timescales up to 1 Gyr. The simulation was
synchronized to two surveys, 1) the Vera C. Rubin Observatory’s Legacy Survey
of Space and Time (LSST) and 2) the Zwicky Transient Facility (ZTF) public and
partnership survey. Upper limits were calculated for the ETA population for
ZTF (with assumed null detection) and LSST’s cadence simulations. For LSST, a
null detection is within three sigma of the simulation results out to H=19 for
the twilight survey and H=20 for the baseline survey. Due to the Yarkovsky
Effect, no ETAs are stable on billion year timescales are likely to be
detected. ETAs large enough to remain stable on billion year timescales are
very rare relative to the rest of the ETA population. A null detection in LSST
will restrict that population to tens of objects larger than 100 meters. The
null detection by ZTF to date has already done so.

Trojan asteroids, Earth Trojans — methods: numerical, N-body Problem

tablenum SIXtablenum

##  1 Introduction

Trojan asteroids111‘Trojan’ is generally used for objects co-orbiting with
Jupiter; however, given the general nature of these minor planets across the
Solar System, we use the term ‘Trojan’ generically and refer to which planet
they co-orbit throughout. are a class of co-orbiting objects that librate
about the L4 or L5 Lagrangian points within a planet’s orbit. Mapping the
distributions of possible Earth Trojan Asteroids (ETAs) is important to guide
observational searches. Librations of ETAs may carry the asteroids nearly the
entire distance from the planet to the L3 point on the far side of the orbit
behind the Sun. In the rotating frame, the gravitational potential lines look
like tadpoles. In this frame, tadpole orbits occupy a half-toroid region,
including either the L4 or L5 point. The region of the sky where ETAs
potentially exist is difficult to observe due to the low Solar elongations.
The geometry also tends to have small reflection angles, making observations
even more difficult from Earth. However, synoptic surveys such as those
carried out by the Zwicky Transient Facility (ZTF, Bellm et al., 2018) and the
Vera C. Rubin Observatory (LSST Science Collaboration et al., 2009; Jones et
al., 2020) could still overcome the difficult observing conditions with
numerous observations at the periphery of the ETA regime. Twilight surveys, in
particular, could be powerful probes into the population. ZTF has been
carrying out a twilight survey for years (see e.g., Ye et al., 2020), and
there is still discussion of whether or not Rubin Observatory’s Legacy Survey
of Space and Time (LSST) will or will not have a dedicated twilight survey
(Seaman et al., 2018).

###  1.1 Currently known ETAs

Currently, there are two known Earth Trojan asteroids (ETAs): 2010 $TK_{7}$
and 2020 $XL_{5}$ (Connors et al., 2011; de la Fuente Marcos & de la Fuente
Marcos, 2021), which both librate far from the Lagrange points and are only
transient in nature and will only be co-orbiting with Earth for a short period
of time compared with the age of the Solar System. Both known ETAs are on
tadpole orbits that librate about Earth’s L4 Lagrange point, traveling much of
the distance between Earth and Earth’s L3 Lagrange point behind the Sun. The
orbits of both ETAs are stable on the order of thousands of years, with 2020
$XL_{5}$ being the more stable of the two. 2020 $XL_{5}$ was shown to likely
have been captured in the 15th century and will remain stable for
approximately four thousand years (Santana-Ros et al., 2022).

###  1.2 Scientific importance of ETAs

There are several potential sourcing mechanisms for Trojan asteroids. First,
they may be formed in the early proto-planetary disk near either a planet’s L4
or L5 Lagrange points (Montesinos et al., 2020). Second, Trojans may be
captured debris produced from major planetary impacts. Finally, they may be
sourced from perturbed populations of asteroids elsewhere in the Solar System,
either during planetary migration or from ejected asteroids from the main-
asteroid belt. The first two mechanisms would have occurred early in the
formation of the Solar System, and we consider these to be “primordial” in
nature. For ETAs, primordial sourcing through the impact hypothesis may also
explain lunar formation. The giant impact hypothesis has been proposed, where
another planet impacted the proto-Earth from a 1 AU orbit in the solar nebula.
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 proto-Earth to form the moon. A cloud of debris from such a
collision may have provided a primordial source of ETAs (Canup, 2004). 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. The
existence of primordial ETAs would offer prime low $\delta v$ targets for
rendezvous missions to study the early inner solar system and especially the
Earth-Moon system (Malhotra, 2019). The NASA mission Lucy mission will soon
begin its 12-year journey to several families of Jupiter Trojans, which will
start to answer some of the scientific questions we have discussed (Zhou et
al., 2020). While the Lucy mission is focused on primordial Jupiter Trojans,
it will shed light on how scientifically rich these objects are, proving that
future Earth Trojan missions will be important in understanding our Solar
System.

There are several potential sourcing mechanisms for Trojan asteroids. First,
they may have settled into a planet’s L4, or L5 Lagrange points directly from
their formation in the proto-planetary disk (Montesinos et al., 2020). Second,
Trojans may be captured debris produced from major proto-planetary impacts in
the inner Solar System. Finally, they may be sourced from perturbed
populations of asteroids elsewhere in the Solar System, either during
planetary migration or from ejected asteroids from the main-asteroid belt. The
first two mechanisms would have occurred early in the formation of the Solar
System, and we consider these to be “primordial” in nature in that their
chemical composition would inform inner Solar System formation processes.
Primordial sourcing through the impact hypothesis may also explain lunar
formation. The giant impact hypothesis has been proposed, where another planet
impacted the proto-Earth from a 1 AU orbit in the solar nebula. 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
proto-Earth to form the moon. A cloud of debris from such a collision may have
provided a primordial source of ETAs, (Canup, 2004). 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. The existence of primordial
ETAs would offer prime low $\delta v$ targets for rendezvous missions to study
the early inner solar system and especially the Earth-Moon system (Malhotra,
2019). The NASA mission Lucy mission will soon begin its 12-year journey to
several families of Jupiter Trojans, which will start to answer some of the
scientific questions we have discussed, (Zhou et al., 2020).

Finally, ETAs have even been posed as prime targets for commercial space
endeavors (Staehle, 1983).

###  1.3 Observational searches for ETAs

Observers have searched for ETAs with a few relatively small astronomical
twilight surveys (Whiteley & Tholen, 1998; Connors et al., 2000; Markwardt et
al., 2020; Lifset et al., 2021). Markwardt et al. (2020); Lifset et al. (2021)
were able to use a null detection to place competitive upper limits on the ETA
population, predicting an L4 ETA population of ${N}_{ET}<1$ for $H=14$,
${N}_{ET}<10$ for $H=16$, and ${N}_{ET}<938$ for $H=22$. Efforts have been
proposed on how to best observe the ETA population utilizing Solar System
exploration missions (Cambioni et al., 2018; Yoshikawa et al., 2018; Malhotra,
2019), ground-based telescopes (Todd et al., 2012, 2014; Seaman et al., 2018;
Ye et al., 2020), and NEO Surveyor Mainzer et al. (2021) which is planned for
a launch this decade to search for potentially hazardous asteroids including
ETAs.

###  1.4 Previous ETA simulations

Numerous studies of ETAs have been carried out in the literature featuring
numerical simulations (Wiegert et al., 2000; Evans & Tabachnik, 2000; Morais &
Morbidelli, 2002; Todd et al., 2012; Lifset et al., 2021; Zhou et al., 2019;
Napier et al., 2022; Yeager & Golovich, 2022a, b). These studies have
typically centered around the questions of sourcing, lifetimes, stability
zones, and observability of the stable population. We refer readers to our
previous paper analyzing the MEGASIM for lifetime and stability of ETAs
including a detailed literature review (Yeager & Golovich, 2022c). Regarding
observability, several authors produced spatial maps of time-averaged ETA
positions to guide observers. For example, Wiegert et al. (2000) created ETA
distribution maps from simulations of two synthetic populations of Earth
Trojans. These models injected Earth co-orbiting asteroids uniformly within
the region between L2 and L3 points. When viewing ETA orbits from Earth, the
resulting peak spatial density is a few degrees displaced toward the Sun from
the L4 and L5 Lagrange points. Evans & Tabachnik (2000) and Morais &
Morbidelli (2002) found stable zones with inclination between 10° and 45° and
with a wide range of ecliptic longitudes encapsulating the L4 and L5 points
but spanning tens of degrees to either side of L4 or L5. Todd et al. (2012)
modeled albedo, size, and orbital distributions in an effort to define an
‘optimal’ search strategy for ETAs. Based on their simulations, the authors
suggest a twilight observing campaign to search the wide swath of area that
the ETAs in their simulation cover. The only known asteroid at the time,
$2010TK_{7}$, would appear in such a survey.

An in-depth study of ETAs that used MEGASIM data can be found in Yeager &
Golovich (2022b). While the currently known ETAs exhibit orbital stability on
the order of thousands of years, objects librating about the Earth’s L4 or L5
Lagrange points can be stable on the order of millions or billions of years.
An ETA’s ability to remain stable depends on several factors, e.g., spin rate,
shape, size, albedo, density, and surface conductivity. These properties play
a role in non-gravitational forces such as the Yarkovsky Effect. Smaller ETAs
experience larger effects from Solar radiation, thus driving them from the co-
orbital regime more quickly. Simulations by (Zhou et al., 2019) find that only
ETAs larger than $\sim$100 meters remain stable for 1 Gyr.

In this paper, we study the MEGASIM data to assess the observability of ETAs
in ongoing and upcoming synoptic surveys. In §2, we give a brief summary of
the simulation details, and we describe our method of determining the flux of
individual asteroids in the simulation. In §3, we present spatial
distributions and discuss translating our simulation output into observational
surveys. In §4, we present the results regarding the observability of ETAs in
two synoptic surveys. Finally, in §5, we offer conclusions of our analysis.

##  2 The MEGASIM

The Multitudinous Earth Greek Asteroid Simulation (MEGASIM), announced first
in Yeager & Golovich (2022a), contains 11.2 million ETA orbits initialized on
trajectories that span beyond the expected ETA stability regime around the L4
Lagrange point. The orbits have now been carried out to a maximum of 4.5
billion years, allowing a study of the longevity of ETA orbits. Here we offer
a summary of those findings, but we refer readers to that work for more
details (Yeager & Golovich, 2022b). We found that the survival of the ETA
population follows a stretched exponential form. Ejections of ETAs were
closely tied to resonances in Earth’s orbit, such as the Milankovitch cycles,
as well as resonances from variations in the orbits of the other planets. This
paper will focus on the spatial distributions and observability of ETAs with
ZTF and LSST.

###  2.1 Initialization and Parallelization

Briefly here we summarize the simulation setup. The MEGASIM data set consists
of two distinct simulations run with different integrators. The work in this
paper focuses on the WHFast data set (Rein & Liu, 2012; Wisdom & Holman, 1991;
Rein & Tamayo, 2015; Rein & Spiegel, 2015). Since the WHFast integrator is
significantly faster, we were able to integrate to the full 4.5 billion year
age of the Solar System. To achieve parallelization, the 11.2 million ETA
trajectories were placed into batches of 500, resulting in 22,400 parallel
Solar Systems that were integrated forward. 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 JD2459145.61625984 (Giorgini et al.,
1996). We utilized the Catalyst, Quartz, and Ruby clusters managed by
Livermore Computing at Lawrence Livermore National Laboratory
222https://hpc.llnl.gov/hardware/platforms/ruby.

Each batch of asteroids were initialized with six orbital parameters selected
from a random distribution, as described below. Each model was then integrated
forward with REBOUND. The integration time step was set to one-day, and
position and velocity information for each simulation was output every 1000
years.

The ETA orbits were initialized with a true longitude ($\theta$) 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 was
randomly assigned values between 0 and $2\pi$, 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
given in Equation 1.

$\boldsymbol{\mu}=\begin{pmatrix}\bar{a}\\\ \bar{e}\\\
\bar{i}\end{pmatrix}=\begin{pmatrix}1\,\text{au}\\\ 0\\\
0\degree{}\end{pmatrix}~{}\boldsymbol{\Sigma}=\begin{pmatrix}0.025\,\text{au}&0&0\\\
0&0.075&0\\\ 0&0&15\degree{}\end{pmatrix}.$ (1)

###  2.2 Modeling ETA intrinsic properties

Every ETA in the MEGASIM is assigned an albedo, asteroid type, diameter, and
absolute magnitude ($H$). The albedo of the ETAs are randomly sampled from the
distribution curve of NEAs determined by Wright et al. (2016) and with a
functional form given by Equation 2

$p(p_{v})={\exp}\left(\dfrac{-p^{2}_{v}}{2d^{2}}\right)\frac{f_{D}p_{v}}{d^{2}}+{\exp}\left(\dfrac{-p^{2}_{v}}{2b^{2}}\right)\,\frac{p_{v}}{b^{2}}\,(1-f_{D})$
(2)

where $p$ is the probability density function, $p_{v}$ is the asteroid albedo,
the dark fraction $f_{D}=0.253$, bright peak $b=0.168$ and dark peak
$d=0.030$. The albedo distribution consists of two Rayleigh distributions with
peaks at a value of 0.030, accounting for 25.3% of the NEAs observed by WISE
and 0.168 for the remaining 74.7% of the NEAs.

The ‘type’ of the asteroid is either an ‘S’ or ‘C’ type and is determined via
a second level of random sampling after the albedo is determined. The two
Rayleigh functions that produce Equation 2 correspond to an ‘S’ or ‘C’ type
asteroid. The height of each function is used to weight the random sample that
determines whether an ETA is of ‘S’ or ‘C’ type. The ETA type determines the
correction magnitude for the final apparent brightness calculation.

The asteroid diameter is obtained from $H$ and albedo ($p_{v}$) using Equation
3.

$D\textsubscript{ast}=\frac{1329000}{\sqrt{p_{v}}}\times 10^{-0.2H}$ (3)

$H$ was randomly sampled from an extended version of the distribution
presented in Granvik et al. (2018), which covered $17<H<25$ with a break in
the power law slope at $H=23$. To sample from the distribution, the number of
asteroids at a given absolute magnitude was approximated using two linear
functions, but we extrapolated the distribution to $H=28$.

The apparent magnitude of an ETA is obtained via Equation 4

$\displaystyle
m_{V}=H+5\log{}_{10}\left(\frac{R_{AS}R_{AE}}{R_{0}^{2}}\right)-2.5\log{}_{10}\left(q(\alpha)\right)$
(4)

where $R_{AS}$ is the distance between the Sun and the asteroid, $R_{AE}$ is
the distance between the asteroid and the center of Earth, $R_{0}$ is a
conversion factor of 1 AU to meters (the unit used for all distances) and
$q(\alpha)$ is the phase integral of a diffuse reflecting sphere Equation 5

$\displaystyle
q(\alpha)=\frac{2}{3}\left(\left(1-\frac{\alpha}{180\degree}\right)\cos{\alpha}+\frac{1}{\pi}\sin{\alpha}\right)$
(5)

where $\alpha$ is the phase angle between the Sun-Asteroid-Earth.

The magnitude of the Sun, $M_{Sun}$ for calculations in this paper is chosen
to be 4.80, the magnitude of the Sun in Johnson V-band (Willmer, 2018).
Corrections from the Johnson V band to both the LSST and ZTF filters are
calculated by DeMeo et al. (2009); Chance & Kurucz (2010); Chesley & Veres
(2017), the values of which are provided in Table 1.

=-0.75in  Class V-u V-g V-r V-i V-z V-y C -1.614 -0.302 0.172 0.291 0.298
0.303 S -1.927 -0.395 0.255 0.455 0.401 0.406

Table 1: Color corrections from Johnson’s V-band to LSST filter system.

##  3 Distributions of ETAs

###  3.1 MEGASIM coordinates

The Cartesian coordinates of the MEGASIM are aligned with aspects of the Solar
System as it was on $t_{0}$=JD2459145.61625984. The X-axis is aligned with the
vernal equinox at simulation time = $t_{0}$, with both the X and Y-axis in the
plane of the ecliptic. The Z-axis is aligned with the ecliptic north pole. We
define a rotating ecliptic spherical coordinate system centered on Earth. The
Sun’s latitude and longitude are always at the origin in this rotating frame,
and the idealized L4 point is at (0°, 60°). This simplifies the analysis by
allowing us to simply stack parallel simulations with one another. It also
satisfies the goal of forecasting the observability of ETAs.

###  3.2 Stacked distribution maps

Figure 1 shows the number density of all ETA positions stacked in 100 kyr
intervals up to 1 Gyr. As ETAs are removed from the MEGASIM due to
instability, they no longer contribute to these maps. Unstable orbits are more
numerous but offer significantly shorter trajectories than the few orbits that
are stable on long timescales. The figure is presented in the ecliptic
rotating coordinate system discussed above. The black circle is the L4 point
(0°,60°), and the red cross is the bin with the highest number of ETAs. The
peak density is 7° Sun-ward offset from L4. The results are largely consistent
with a similar figure from Wiegert et al. (2000) but with a much higher
resolution. The stacked 2d histograms reveal the large on-sky area over which
ETAs can travel, but for the most part, ETAs are found within 10° of the
ecliptic. There are several bands of higher and lower stability. These are
evident in later figures with different coordinates presented. Notable are
ridges and gaps above and below the ecliptic plane as longitudes between 20°
and 50°, as well as wings that flare away from the ecliptic at longitudes
beyond L4.

Figure 1: Stacked map of ETAs from the WHFast integrated simulation of all
ETAs every 100 kyr up to 1 Gyr. The ecliptic coordinates are explained in
Section 3.1. The Sun is at the origin, the black dot is L4 at (0,60°) and the
red cross indicates the location of highest density, which is 7° offset from
L4. The bin amounts are scaled to the bin with the highest density.

Weighting the stacked ETA trajectories of Figure 1 by their apparent
brightness has only a minor effect in determining where they may be observed.
If the ETA on-sky positions are weighted by their flux, the location of the
peak density is only shifted by about 0.5°. One dimensional cuts of the
ecliptic longitude and latitude of Figure 1 are shown in Figure 3. The notable
features are a steep fall off as the ecliptic latitude increases and a local
minimum at 40° and local maximum around 35° ecliptic longitude. In Figure 2 we
show the difference obtained from the ETA trajectories by brightness, binning,
and normalizing all bins such that the heights exist in the range of 0 to 1,
we then subtract off the bin heights of Figure 1. ETAs within 10° of the
ecliptic interior to L4 are generally fainter than the density map would
suggest (hence the blue colors interior to the L4 point). This is due to the
majority of ETAs in this portion of the sky existing at fairly far distances
from the Earth, resulting in lower apparent brightness. Those found in the
region exterior to L4 from the Sun tend to be close to their nearest approach
to Earth on their trajectory. The phase angle in the brightest region does
result in a reduction of apparent magnitude, but the effect is sub-dominate to
the proximity effect with Earth. For observing considerations, the region
behind L4 is much easier to observe due to the greater solar elongation
compared to trying to observe the generally fainter ETAs interior on the sky
to L4.

Figure 2: A difference map produced by weighting ETA trajectories by apparent
magnitude, binning into a 2d histogram, normalizing all bins such that the
heights exist in the range of 0 to 1 then subtracting off the bin heights of
Figure 1. The yellow colors indicate where the ETA apparent magnitude
weighting is positive, green is near zero, and purple is negative. The black
circle is the L4 point, and the red cross shows the location of the highest
deviation from the density plot in Figure 1. The peak of the flux-weighted
density map is located 6.5° Sun-ward from L4. Figure 3: One degree slices of
figure 1 are made in both latitude and longitude. The y-axis is the percent of
ETAs trajectories found for a given one degree thick cross section in figure
1.

Figure 4 shows the projected density of stacked ETAs positions in the ecliptic
plane. It reveals that ETA trajectories extend into the Solar System as far as
the orbit of Venus and span the full range of longitudes from Earth to L3. As
one would expect, the highest ETA densities are found around L4, but slightly
counterclockwise from L4, which agrees with the location indicated in Figure
1. Measuring an angle counter-clockwise from Earth at the top of Figure 4, the
highest density is located about 5 degrees further from Earth. ETA density
around the peak remains within 10% of the peak as close as 45° out to about
80° from Earth. From 90° to 105°, there is a region of reduced ETA density
before another local maximum in ETA density found near 115-120° from Earth.

Figure 4: A top-down view of the Solar System. The Sun is at (0,0), and the
concentric black circles starting nearest the Sun indicate the approximate
orbits of Mercury, Venus, Earth, and Mars. Radial dashed lines are provided
for angular reference at 30°, 45°, 60°, 75°, 90°, 105°, 120°, 135°, 150°
measured counter-clockwise from the location of Earth in the figure. Binning
is the same process as done for Figure 1

###  3.3 Distribution of billion year-stable ETAs

Figure 5 shows single snapshots in time and only includes ETAs that have
survived at least 1 Gyr of the MEGASIM integration. The top panels are the
initial positions of all billion year (Gyr) stable ETAs, and the bottom panels
are 100 Myr later. To state the most obvious first, Gyr-stable ETAs are of low
inclination and do not extend beyond 10°ecliptic latitude. The spatial
distribution, unsurprisingly, does not span far from that of Earth’s orbit.
The same density ridges and gaps seen in Figure 4 persist in the Gyr stable
orbits.

The tail of the ETA distribution is within the drawn circular Earth orbit at
100 Myr, seen in the bottom right panel of Figure 5. Other snapshots in time
reveal the distribution of ETAs oscillates in and outside of the drawn
circular orbit of Earth at a frequency higher than the simulation output of
1000 years. This indicates that there is a bulk forcing of the ETA population
throughout their libration trajectories.

Figure 5: The distributions of only those ETAs that survive at least 1 Gyr of
simulation. In each panel, the black dot is the approximate location of L4,
the red cross is the bin with the highest Gyr ETA density, and the circles are
at the semi-major axis distances of Mercury, Venus, Earth, and Mars. The
location of the highest Gyr ETA density oscillates between 45°-70° latitude
on-sky and in the X-Y plane, travels as close as 45° leading Earth to as far
as 90°.

###  3.4 The LSST and ZTF surveys

The survey data used for the analysis in the section 4 is the LSST v2.2
baseline ‘lsst_no_twilight_neo_v2.2’ and twilight
‘twi_neo_repeat4_riz_v2.2_10yrs’ OpSim
runs333https://www.lsst.org/scientists/simulations/opsim as well as all ZTF
observations spanning dates 2018-03-17 to 2022-06-01, (LSST project, 2022;
Bellm et al., 2018). The ZTF data set was filtered down to only include
30-second exposures taken by the Palomar P48 observing system, including both
public and private partnership observations. The data set consists of 706,459
30-second observations. The number of observations in the LSST baseline OpSim
was 2,259,570 and twilight 2,078,065.

###  3.5 Translating MEGASIM orbits to survey pointings

Of the 11.2 million MEGASIM initialized orbits, only those at least stable for
15 thousand years were propagated through the survey simulations. This
resulted in the propagation of 368,287 ETA orbits with a maximum time step of
30 seconds and integrated to the exact time of each LSST and ZTF observation.

To simulate how many MEGASIM ETAs could be detected in the LSST and ZTF
surveys, we propagated the initialization simulation from $t_{0}$ to the start
date of the respective surveys. For LSST, forward propagation was done to sync
with the start date of the LSST OpSims, and backward propagation was done to
sync the Solar System with the beginning of the ZTF survey. The shifting of
initialization resulted in on-sky errors of the order $10^{-6}$° in right
ascension (RA) and
$6\text{\times}{10}^{-7}\text{\,}\mathrm{\SIUnitSymbolDegree}$ in declination
(Dec.). Over the entire course of integration, a maximum error of 0.25° in
R.A. and 0.07° in Dec. was reached. Considering the effective field of view
diameter of 1.75° for LSST and 3.4° for ZTF, the deviations between the
simulated and true positions do not significantly affect results.

Integration of the ETA orbits is carried out with a 30-second time step,
except in cases where observation times are shorter than the time step. The
simulation is carried out to the exact timestamp of all observations for the
LSST and ZTF surveys. Once the MEGASIM Solar System reaches the timestamp of a
given survey observation, a check is done to determine which ETAs are at an
R.A. and Dec. within the observations field-of-view. To simplify this check, a
circle is used to approximate the field-of-view of LSST and a square for ZTF.
The effective circle for LSST is taken to be of radius 1.75°, and for ZTF, the
square has side lengths of 3.428° and an area of 47 sq degrees.

ETAs found within a survey’s observation are then cataloged. A further check
is then done to determine if the ETA in the field of view is bright enough to
be detected by the telescope. If the ETA’s apparent magnitude is brighter than
the 5$\sigma$ sensitivity for that exposure, it is cataloged as a detection.

###  3.6 Generating detection statistics

To place statistical weights on ETA detections in the LSST and ZTF surveys.
Detections in each survey simulation were recalculated by assigning all ETAs a
new $H$, diameter, and albedo. The new ETA properties are obtained via
resampling from the $H$ and $p_{v}$ distributions laid out in section 2.2.

This allowed ETA resampling and detections in post without needing to re-
simulate all of the ETA trajectories. The resampling was completed 1000 times
for each LSST Opsim and the ZTF pointings. In each step, a unique detection
list was created following the procedure described above.

##  4 Results: The detection of ETAs in ZTF and LSST surveys

###  4.1 Detection statistics and cumulative counts

The r-band filter for LSST and the g-band filter for ZTF are the most
successful at detecting ETAs; see Figure 6. ZTF’s g-band is detects more ETAs
than LSST’s g-band. Of the 706,459 ZTF pointings used (38.57%, 272,541) were
in g-band, (56.99%, 402,632) r-band, and (4.43%, 31,286) i-band. For the LSST
baseline OpSim, the distribution of filters used (u, g, r, i, z, y) are:
(18.67%, 387,928), (21.71%, 451,222), (22.66%, 470,990), (19.61%, 407,434),
(10.96%, 227,705), (6.39%, 132,786), respectively. The ratio of the number of
g-band and r-band observations by LSST is of unity but provides an order of
magnitude fewer detections than ZTF.

Figure 6: Detected ETAs by survey filter type. Bins are normalized to the
total propagated 1 Myr stable ETA trajectories (368,287).

The cumulative number of ETA detections for ZTF and the LSST OpSims is
provided in Figure 7. The rate ZTF detects ETAs is roughly the same as the
baseline LSST OpSim. Twilight LSST captures approximately four times the
number of ETA detections as either LSST baseline or ZTF. Figure 7

In Figure 8, we show the relative detection rate for each survey by month.
Variations occur in the detection rate of ETAs as the time a telescope can be
on target toward L4 varies over the year Whiteley & Tholen (1998). Beyond the
variable detection rate, ZTF finds most ETAs from approximately June to
December, whereas LSST makes ETA detections from January to June. These
differences are due to the different hemispheres of each observatory, which
changes the time of year L4 is visible. LSST detects the most ETAs in the
month of April. ZTF detects nearly the same number of ETAs from August to
November.

Figure 7: The cumulative # of unique ETA detections made by each survey for
the fiducial simulation. Figure 8: The normalized number of detections by
month made by each survey. Figure 9: The cumulative % of unique ETA detections
made by each survey for the fiducial simulation. The total population of
unique ETAs is 368,287.

The number of distinct detections per asteroid for each survey is shown in
Figure 10. There is a stark difference between the two LSST surveys. This is
because the twilight survey observes lower Solar elongations nightly, which is
capable of recovering the same asteroid thousands of times as it slowly
librates through the twilight survey footprint. The ZTF survey also has a
twilight survey, which explains the few instances of an ETA being detected
thousands of times; however, it is far less sensitive, so this only works for
a few bright ETAs. The baseline LSST survey recovers many asteroids numerous
times, but not nearly as many as the LSST Opsim with a twilight survey. Note
that even the baseline LSST survey still recovers asteroids as many as 2000
times. This is far more than typical near-Earth asteroids are recovered. Jones
et al. (2018) estimated that typical NEOs are recovered $\sim$100 times. The
discrepancy here is that ETAs move through the survey footprint over years
rather than weeks.

Figure 10: Asteroid recovery by each survey.

Figure 11 shows the distribution of $H$ for detected ETAs. Unique detections
by ZTF are few and do not occur for $H>20.5$. For LSST, there is very little
difference in the distribution of $H$ when comparing the baseline and twilight
OpSim surveys. LSST makes the bulk of detections in the range $21<H<25$,
though some detections are made up to $H<28$, which is the limit of the
simulation. This suggests that ETAs come close enough to Earth to be detected
by LSST to extremely faint magnitudes, which is a statement of the Rubin
Observatory’s sensitivity more than anything.

Figure 11: Detected ETA % by $H$. Bins are normalized to the total propagated
1 Myr stable ETA trajectories (368,287).

The apparent magnitudes of the ETAs detected by LSST baseline (irrespective of
filter), twilight and ZTF are provided in Figure 12. ZTF cannot detect any
ETAs beyond an apparent magnitude of 21, while most LSST observations are at
apparent magnitudes between 21 and 23. Detections of ETAs quickly fall off
after an apparent magnitude of 23 but are still possible by LSST out to 24.3.
The LSST twilight survey is unsurprisingly much more successful at repeatedly
detecting ETAs at $m_{app}>21$. The LSST baseline and twilight are similar in
detections of unique ETAs, but twilight produces more than an order of
magnitude more repeat detections.

Figure 12: Total detections of ETAs by apparent magnitude.

Figure 13 shows the proper motions of detected ETAs peak at 0.05 ”/s with an
extremely sharp fall-off tail reaching 0.2 ”/s. The LSST baseline made the
fastest ETA detection of 0.3 ”/s. Though, the twilight survey’s overall proper
motion distribution is very similar to that of the baseline. The on-sky speed
is typical of ETAs at their closest approach to Earth. ZTF, with a larger
field of view, still failed to capture any ETAs moving faster than 0.10 ”/s.
This speaks more to the apparent magnitude limit of ZTF than LSST; if ZTF
could make more detections, a similar small fraction would be expected to
proper motions of 0.3 ”/s.

Figure 13: Total detections ETAs proper motions by the survey.

Figure 14 shows the distance an ETA is from Earth at the time of detection.
ZTF can make ETA detections out to 1.7 AU, the practical limit due to the
angle of solar elongation becoming too small for further distances. As seen in
the X-Y map of ETA trajectories in Figure 4, ETAs travel out to L3; behind the
Sun as viewed from Earth, though detections cannot be made from Earth beyond a
distance of 1.75 au as the line of sight is too near the Sun and neither LSST
nor ZTF point at solar elongations less than 36°. ZTF detects 5-10x fewer ETAs
at all distances compared to LSST’s baseline and twilight survey. The LSST
baseline is nearly 1:1 with the twilight survey for ETA detections of 0.75 au
and closer. However, a significant deviation in detections occurs after 0.75
au, where the LSST twilight survey becomes much more successful at detecting
ETAs. The ZTF survey probes a few degrees closer to the Sun than the LSST
twilight survey allowing a few more detects at distances out to 1.75 au.

Figure 14: Total detections of ETAs by distance from the Earth in units of au.

Figure 15 shows three different views (one in each row) of the ETA detection
distributions for the three surveys (one in each column). The top row shows
the distribution of ETA detections viewed from the pole of the Solar System.
The middle row shows the distribution as viewed from Earth, where (0,0) is the
position of the Sun and (0,60) is the L4 point. The bottom row shows the R.A.
and Dec. of the ETA detections. All three surveys can detect ETAs within
similar X-Y bounds, as seen in the top row of figures. The detections made by
ZTF all occur closer to the orbit of Earth, with approximately an even number
of observations made over all distances from Earth. This is due to ZTF only
picking up a few relatively bright ETAs but observing them repeatedly over a
large libration about L4. Comparing the LSST baseline to the twilight results,
a stark difference is seen in the X-Y detection distribution as the twilight
survey successfully detects many more ETAs at solar elongations of 30-35
degrees, which correspond to distances from Earth of 0.75 au to 1.7 au. The
increased bright region in Figure 15 middle figures of the first two rows
corresponds to the bump of twilight detections in Figure 14. The ability of a
survey to make more observations at solar elongations within L4 significantly
increases detections of ETAs. The R.A. and Dec. distribution follows the
ecliptic but also re-highlights that ZTF observes ETAs between June and
December, and LSST, because of its southern latitude, observes the L4 ETAs
from January to late June.

Figure 15: 2d histograms of detected ETA distributions. Binning is 50x50. The
three columns from left to right are for ZTF, LSST twilight, and LSST baseline
surveys. The three rows are: 1) X vs. Y; top-down view of the solar system
with the Sun at (0,0) and from the left in the approximate orbits of Mars,
Earth, Venus, and Mercury. The frame is such that the position of Earth is
fixed to the point (-1,0). The markers for the other planets are not their
true position but denote the orbits. 2) Latitude vs. longitude; this is an on-
sky coordinate system viewed from Earth, the same as used in Figure 1. The Sun
is fixed at (0°,0°), and L4 is at (0°,60°), denoted by a red dot. 3) R.A. vs.
Dec.

###  4.2 The detectable fraction of the ETA population by LSST and ZTF

The fraction of ETAs detected was determined for 1000 bootstrap samples of the
ETA populations for each survey. In Figure 16, the survey completeness for ETA
detections as a function of $H$ is shown. The gray line represents the median
simulated ETA $H$ distribution, with 68% confidence interval (CI) shown as
error bars. The median survey completeness is shown as a solid color lines.
The 68% confidence intervals are shown as dotted lines. Low number statistics
at the bright end of the $H$ distribution tend to cause an all-or-nothing
detection of ETAs, giving the wide zero to one error bars for the survey
completeness. The top left panel shows that the median number of detected ETAs
for each survey is zero at $H<18$. This is due to the fact that brighter ETAs
are very rare but does not suggest that LSST and ZTF would fail to detect such
bright ETAs but rather that the occurrence of a bright ETA within the survey
footprint is rare in our bootstrap sampling procedure. As the simulation time
progresses, the low number statistics become more prominent for smaller and
smaller (originally more numerous) ETAs. In the bottom right panel, no ETAs
with $H<19$ are present in any of the bootstrapped populations.

The subpanels of Figure 16 show the fractional completeness for a series of
different stability subgroups of ETAs that were selected by their lifetime in
the MEGASIM. ETAs were removed from the MEGASIM if they crossed the plane
perpendicular to the ecliptic that intersects the Earth and L3. In the
MEGASIM, only gravitational forces were used to compute the motion of the
ETAs. Additional forces, such as thermal effects like the Yarkovsky effect,
would, in reality, further reduce the number of ETA trajectories that remain
on long timescales especially for smaller asteroids. For orbits that persist
on Myr timescales, additional forces are negligible. If an orbit persists on
the order of 1 Gyr, ETAs with smaller diameters ($\sim 100$ meters) will be
pushed out of tadpole stability regime Zhou et al. (2019). The corresponding
absolute magnitude for 140 meter asteroids is approximately $H=20.1$ for
S-type and $H=22.2$ for C-type ETAs. In the bottom right panel, the only
detections that remain are fainter than these limits.

Figure 16: The fraction of simulated ETAs detected in each survey. In green is
the LSST v2.2 baseline, orange is the LSST twilight, and blue is ZTF. The gray
line represents the median ETA absolute magnitude distribution, with a 68%
confidence interval (CI) shown as error bars. The fractional survey
completeness of 68% CI is shown as dotted lines.

###  4.3 Estimating the upper limit of the ETA population

Figure 17 is the cumulative upper limit for the ETA population determined via
the fractional completeness curve. Following the method of Markwardt et al.
(2020), we determine the upper limit population allowable assuming a null
detection of ETAs in each survey. We calculate the upper limit using Equation
6, where $U(H)$ is the upper limit at a given absolute magnitude, $f(H)$ is
the fraction of ETAs detected at a given absolute magnitude, and the 3 in the
numerator implies the upper limit is within three standard deviations of the
null detection.

$U(H)=\frac{3}{f(H)}$ (6)

Figure 16 shows $U\left(H\right)$ for the same four stability regimes shown
previously. In Figure 16, there are portions of the completeness curves that
are zero, which results in infinite error bars. To handle this, undefined
lower and/or upper error bars caused by division by zero are replaced by -3 or
+3, respectively.

###  4.4 The Yarkovsky Effect and ETA Detections

In the MEGASIM, only gravitational forces were used to propagate the ETA
orbits. Here we estimate the number of gravitationally stable ETAs that may be
driven unstable by non-gravitational effects. The dominant additional force to
consider is the Yarkovsky Effect (e.g., Chesley et al., 2003), which causes
asteroids to rotate and increase their semi-major axis due to Solar radiation
pressure. Simulations of ETA orbits with the Yarkovsky Effect included were
done by Zhou et al. (2019). The Yarkovsky Effect’s strength depends on the
shape, albedo distribution, surface conductivity, spin, and density of a given
asteroid. The additional force from the Yarkovsky Effect was found to remove
ETAs within 1 Gyr for ETAs with diameters smaller than 130 m and 90 m for
prograde and retrograde rotation, respectively (Zhou et al., 2019).

To estimate the changes to our survey completeness curves that the Yarkovsky
Effect would impose, a cut-off on the simulated ETA population was made using
both the lifetime and diameter of a given ETA. To estimate the removal rate of
gravitationally stable ETAs in the MEGASIM, we use Equation 7.

$D_{unstable}=\frac{lifetime\,\,\,100\,\,meter}{1\,\,Gyr}$ (7)

where the lifetime is the time an ETA remained bound to L4 in the MEGASIM in
Gyr. An ETA is assumed to be removed from the population and, therefore, not
detectable if the diameter of an ETA is less than ($D_{unstable}$). The
diameter of 100 m was chosen as a round value between the diameter estimated
by Zhou et al. (2019) for ETAs with prograde and retrograde spins, which are
removed by 1 Gyr.

The results of these cuts on the fraction of the population detected are shown
in Figure 18. The x-axis is the absolute magnitude ($H$) of the detected ETAs,
and the y-axis is the fraction of the seemingly detected ETAs in Figure 16
that would be removed via the Yarkovsky Effect at given times. The Yarkovsky
Effect is not strong enough to drastically change our estimated survey
detections for ETAs stable for less than 100 million years. All detected ETAs
with lifetimes greater than 1 Gyr are no longer possible if the Yarkovsky
Effect is considered. ZTF is not included because there are no timescales
where ZTF detections remain after this procedure. We reiterate here that this
is due to small number statistics in our bootstrap populations in the
detection efficiency analysis above.

##  5 Conclusion

Figure 17: Cumulative upper limit on ETA populations for the three surveys.
The dotted lines show the 68% CI. Figure 18: The percent of detected ETAs with
lifetimes greater than a given age lost per absolute magnitude bin. The solid
and dashed dot lines are median values, and the dashed lines are the 68% CI of
the 1000 runs. From top to bottom are for ETAs of ages greater than 100, 250,
and 500 Myr. All detected ETAs with lifetimes of at least 750 Myr are removed.
ZTF is not plotted, as no ETAs were detected once the Yarkovsky effect was
considered.

The work presented here can be categorized into three areas. 1) The spatial
distributions of MEGASIM ETA trajectories, 2) the detectability of MEGASIM
ETAs by LSST and ZTF, and 3) an upper limit of an actual ETA population
assuming a null detection.

###  5.1 ETA spatial distributions

When viewing L4 ETAs from Earth, the highest on-sky concentration is found
6.5° closer to the Sun than the L4 Lagrange point. 30% of ETA trajectories are
found within 45-60° solar elongation and another 30% are found beyond 60°. 70%
of all ETAs exist within 5° of the ecliptic plane across all longitudes.
Surveys of ETAs from the ground are challenging inside L4. Detections cannot
be made at all times of the year due to the Earth’s axial tilt (Whiteley &
Tholen, 1998), and observations require longer exposure times due to higher
sky brightness at twilight and higher airmass than typical observations at
opposition.

The distribution of ETAs as viewed looking down the pole of the ecliptic plane
shows a toroidal volume that tapers near the Earth and L3 Lagrange point. ETAs
can librate in as far as the orbit of Venus and out to 1.3 au especially if
they are only stable for short periods of time. Longer stable ETAs mostly stay
within 0.98 - 1.02 au of the Sun. In the ecliptic plane, the highest density
of ETAs oscillates over a considerable distance. High densities occur along a
ridge 45° leading Earth’s orbit to 90° with the peak within a few degrees of
the L4 point if averaged over time. The sloshing of the entire population
takes long periods of time, which raises the possibility that modern surveys
for the population are unlucky with the population anomalously far from L4 at
the present day.

As the lifetime of ETAs increases, the region that surviving ETAs travels over
tightens onto the Earth’s orbit. Long-lived ETAs do not approach Earth or L3
as closely and remain in a tighter bunch around L4. This is one reason for
fewer detections for longer-term stable ETAs in Figure 16, even for ZTF and
the LSST twilight survey, which are notionally designed for such low Solar
elongation observations.

The on-sky location where ETAs are brightest does not correlate to the most
likely position on the sky. The brightest on-sky areas are found at solar
elongations outside of L4 at angles of 60° to 75°. The brightest location is
on the ecliptic, but elevated brightness regions spread beyond L4 to latitudes
of $\pm 20$. This is apparent in Figure 2 as a flaring yellow regime to the
right of L4, and it is due to the improved observing geometry.

###  5.2 Ability of LSST and ZTF to detect ETAs

Due to the latitude that ZTF and LSST are located, the time of year each can
detect L4 ETAs differs (see Figure 8). ZTF finds most ETAs from June to
December and LSST from January to June. In our simulations, ZTF was
unsuccessful at detecting the vast majority of the ETA population due to
limited photometric sensitivity ($H<19$ for most detections). Of the ZTF
detections, only 18 unique ETAs were observed. Though, because of the field of
view of the ZTF telescope repeat detections are on the order of thousands,
which indicates that any detected ETA by ZTF would be very well characterized
and tracked.

Both the LSST baseline and twilight survey recovered similar numbers of unique
ETAs (332 and 342, respectively). The twilight survey, however, was able to
make several thousand repeat detections, with three of the asteroids observed
over 10,000 times. These results show that a twilight survey will not increase
the likelihood of individual ETA detections drastically compared to the
baseline survey. Though, a twilight survey would provide us with several
orders of magnitude more observations of select ETAs. Observing asteroids on
the order of thousands to tens of thousands of times is unprecedented and
would be useful for constraining characterizing individual asteroids.

###  5.3 The possible ETA population

ZTF did not have a statistically significant number of detections to provide
meaningful constraints on the upper limit of the ETA population. However, both
the LSST baseline and twilight surveys provided enough detections to estimate
an upper limit for ETAs of Myr up to Gyr stability times. Due to the Yarkovsky
Effect, C-type and S-type ETAs of absolute magnitudes greater than $H=22$ or
$H=24$ will be driven unstable within 1 Gyr. The rarity of ETAs of sizes
greater than 100 m in diameter paired with the fact that smaller asteroids
will be driven unstable, our results indicate that a null detection by LSST
will mean that there are no remaining Gyr-stable or primordial ETA
populations. On the other hand, the existing ETA discoveries are indicative of
a small and transient population. A null detection in LSST will restrict that
population to tens of objects larger than 100 meters. The null detection by
ZTF to date has already done so.

We would like to thank Bryce Bolin for compiling the ZTF pointings information
and helping us decipher the catalog for our purposes here. We thank Michael
Schneider for valuable discussions. This work was performed under the auspices
of the U.S. Department of Energy by Lawrence Livermore National Laboratory
under Contract DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program
under Project 20-ER-025. Computing support for this work came from the
Lawrence Livermore National Laboratory Institutional Computing Grand Challenge
program. This material is based upon work supported in part by the National
Science Foundation through Cooperative Agreement AST-1258333 and Cooperative
Support Agreement AST-1202910 managed by the Association of Universities for
Research in Astronomy (AURA), and the Department of Energy under Contract No.
DE-AC02-76SF00515 with the SLAC National Accelerator Laboratory managed by
Stanford University. Additional Rubin Observatory funding comes from private
donations, grants to universities, and in-kind support from LSSTC
Institutional Members. This work was based on observations obtained with the
Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar
Observatory as part of the Zwicky Transient Facility project. ZTF-II is
supported by the National Science Foundation under Grant No. AST-2034437 and a
collaboration including Caltech, IPAC, the Weizmann Institute for Science, the
Oskar Klein Center at Stockholm University, the University of Maryland/Joint
Space-Sciences Institute, Deutsches Elektronen- Synchrotron and Humboldt
University, Lawrence Livermore National Laboratory, TANGO/NTHU, Taiwan, the
University of Wisconsin at Milwaukee, Trinity College, Dublin, and IN2P3,
France. Operations are conducted by COO, IPAC, and the University of
Washington.

## References

  * Belbruno & Gott (2005) Belbruno, E., & Gott, J. Richard, I. 2005, AJ, 129, 1724, doi: 10.1086/427539
  * Bellm et al. (2018) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018002
  * Cambioni et al. (2018) Cambioni, S., Malhotra, R., Hergenrother, C. W., et al. 2018, in Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 1149
  * Canup (2004) Canup, R. M. 2004, Icarus, 168, 433, doi: https://doi.org/10.1016/j.icarus.2003.09.028
  * Chance & Kurucz (2010) Chance, K., & Kurucz, R. L. 2010, Journal of quantitative spectroscopy and radiative transfer, 111, 1289
  * Chesley & Veres (2017) Chesley, S. R., & Veres, P. 2017, arXiv preprint arXiv:1705.06209
  * Chesley et al. (2003) Chesley, S. R., Ostro, S. J., Vokrouhlický, D., et al. 2003, Science, 302, 1739, doi: 10.1126/science.1091452
  * Connors et al. (2000) Connors, M., Veillet, C., Wiegert, P., Innanen, K., & Mikkola, S. 2000, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 32, AAS/Division for Planetary Sciences Meeting Abstracts #32, 14.07
  * Connors et al. (2011) Connors, M., Wiegert, P., & Veillet, C. 2011, Nature, 475, 481, doi: 10.1038/nature10233
  * de la Fuente Marcos & de la Fuente Marcos (2021) de la Fuente Marcos, C., & de la Fuente Marcos, R. 2021, Research Notes of the American Astronomical Society, 5, 29, doi: 10.3847/2515-5172/abe6ad
  * DeMeo et al. (2009) DeMeo, F. E., Binzel, R. P., Slivan, S. M., & Bus, S. J. 2009, Icarus, 202, 160, doi: https://doi.org/10.1016/j.icarus.2009.02.005
  * Evans & Tabachnik (2000) Evans, N. W., & Tabachnik, S. A. 2000, Monthly Notices of the Royal Astronomical Society, 319, 80, doi: 10.1046/j.1365-8711.2000.03761.x
  * Giorgini et al. (1996) Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al. 1996, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 28, AAS/Division for Planetary Sciences Meeting Abstracts #28, 25.04
  * Granvik et al. (2018) Granvik, M., Morbidelli, A., Jedicke, R., et al. 2018, Icarus, 312, 181, doi: https://doi.org/10.1016/j.icarus.2018.04.018
  * Jones et al. (2018) Jones, R. L., Slater, C. T., Moeyens, J., et al. 2018, Icarus, 303, 181, doi: 10.1016/j.icarus.2017.11.033
  * Jones et al. (2020) Jones, R. L., Bannister, M. T., Bolin, B. T., et al. 2020, arXiv preprint arXiv:2009.07653
  * Lifset et al. (2021) Lifset, N., Golovich, N., Green, E., Armstrong, R., & Yeager, T. 2021, AJ, 161, 282, doi: 10.3847/1538-3881/abf7af
  * LSST project (2022) LSST project. 2022, Opsim Run Select
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201.  https://arxiv.org/abs/0912.0201
  * Mainzer et al. (2021) Mainzer, A., Abell, P., Bauer, J., et al. 2021, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 53, AAS/Division for Planetary Sciences Meeting Abstracts, 306.16
  * Malhotra (2019) Malhotra, R. 2019, Nature Astronomy, 3, 193, doi: 10.1038/s41550-019-0697-z
  * Markwardt et al. (2020) Markwardt, L., Gerdes, D. W., Malhotra, R., et al. 2020, MNRAS, 492, 6105, doi: 10.1093/mnras/staa232
  * Montesinos et al. (2020) Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, A&A, 642, A224, doi: 10.1051/0004-6361/202038758
  * Morais & Morbidelli (2002) Morais, M., & Morbidelli, A. 2002, Icarus, 160, 1, doi: https://doi.org/10.1006/icar.2002.6937
  * Napier et al. (2022) Napier, K. J., Markwardt, L., Adams, F. C., Gerdes, D. W., & Wen Lin, H. 2022, Planetary Science Journal, 3, 121, doi: 10.3847/PSJ/ac6958
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424, doi: 10.1093/mnras/stu2164
  * Rein & Tamayo (2015) Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376, doi: 10.1093/mnras/stv1257
  * Santana-Ros et al. (2022) Santana-Ros, T., Micheli, M., Faggioli, L., et al. 2022, Nature Communications, 13, 447, doi: 10.1038/s41467-022-27988-4
  * Seaman et al. (2018) Seaman, R., Abell, P., Christensen, E., et al. 2018, arXiv e-prints, arXiv:1812.00466.  https://arxiv.org/abs/1812.00466
  * Staehle (1983) Staehle, R. L. 1983, Astronautics Aeronautics, 21, 44
  * Todd et al. (2012) Todd, M., Tanga, P., Coward, D. M., & Zadnik, M. G. 2012, MNRAS, 420, L28, doi: 10.1111/j.1745-3933.2011.01186.x
  * Todd et al. (2014) —. 2014, MNRAS, 437, 4019, doi: 10.1093/mnras/stt2223
  * Whiteley & Tholen (1998) Whiteley, R. J., & Tholen, D. J. 1998, Icarus, 136, 154, doi: 10.1006/icar.1998.5995
  * Whiteley & Tholen (1998) Whiteley, R. J., & Tholen, D. J. 1998, Icarus, 136, 154, doi: https://doi.org/10.1006/icar.1998.5995
  * Wiegert et al. (2000) Wiegert, P., Innanen, K., & Mikkola, S. 2000, Icarus, 145, 33, doi: 10.1006/icar.2000.6339
  * Willmer (2018) Willmer, C. N. A. 2018, The Astrophysical Journal Supplement Series, 236, 47, doi: 10.3847/1538-4365/aabfdf
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, The Astronomical Journal, 102, 1528, doi: 10.1086/115978
  * Wright et al. (2016) Wright, E. L., Mainzer, A., Masiero, J., Grav, T., & Bauer, J. 2016, The Astronomical Journal, 152, 79, doi: 10.3847/0004-6256/152/4/79
  * Ye et al. (2020) Ye, Q., Masci, F. J., Ip, W.-H., et al. 2020, AJ, 159, 70, doi: 10.3847/1538-3881/ab629c
  * Yeager & Golovich (2022a) Yeager, T., & Golovich, N. 2022a, Research Notes of the AAS, 6, 68, doi: 10.3847/2515-5172/ac62da
  * Yeager & Golovich (2022b) —. 2022b, arXiv preprint arXiv:2206.01837
  * Yeager & Golovich (2022c) —. 2022c, MEGASIM: Lifetimes and Resonances of Earth Trojan Asteroids – The Death of Primordial ETAs?, arXiv, doi: 10.48550/ARXIV.2206.01837
  * Yoshikawa et al. (2018) Yoshikawa, M., Tsuda, Y., Watanabe, S., et al. 2018, in Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 1771
  * Zhou et al. (2019) Zhou, L., Xu, Y.-B., Zhou, L.-Y., Dvorak, R., & Li, J. 2019, A&A, 622, A97, doi: 10.1051/0004-6361/201834026
  * Zhou et al. (2020) Zhou, L., Zhou, L.-Y., Dvorak, R., & Li, J. 2020, A&A, 633, A153, doi: 10.1051/0004-6361/201936332