MEGASIM: Distribution and Detection of Earth Trojan Asteroids

Using N-body simulation results from the MEGASIM data set, 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 start times and timelines of two observational astronomy surveys: (1) the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST) and (2) the Zwicky Transient Facility (ZTF). We calculate upper limits for the number of ETAs potentially observable with both the ZTF and LSST surveys. Due to the Yarkovsky Effect, we find no stable ETAs on billion-year timescales likely to be detected by any ETA survey, as no C-type or S-type ETAs (with H < 22 and H < 24, respectively) are likely to be stable on billion-year timescales, and ETAs large enough to remain stable on billion-year timescales are very rare relative to the rest of the ETA population. We find that a twilight ETA survey will not drastically increase the likelihood of individual ETA detection, but it would provide orders of magnitude more observations of select ETA populations. The null detection to date from ZTF restricts the potential ETA population to hundreds of objects larger than 100 m (at H ≈ 22), while a null detection by LSST will further restrict the ETA population to tens of objects larger than 100 m.


INTRODUCTION
Trojan asteroids1 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).

Currently known ETAs
Currently, there are two known Earth Trojan asteroids (ETAs): 2010 T K 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).

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 δ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 for-mation 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 δ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).

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(Todd et al. , 2014;;Seaman et al. 2018;Ye et al. 2020), andNEO Surveyor Mainzer et al. (2021) which is planned for a launch this decade to search for potentially hazardous asteroids including ETAs.

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°a nd 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, 2010T K 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 ∼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.

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.

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 Laboratory2 .
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 (θ) 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π, assuming a uniform distribution.The semi-major axis (a), ec-centricity (e > 0), and inclination (i) were sampled from a three-dimensional Gaussian with mean and covariance given in Equation ( 1). (1)

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) 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).
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) 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(α) is the phase integral of a diffuse reflecting sphere Equation ( 5) where α 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.

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.

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°S un-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.
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 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.

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.

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 runs3 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 30second 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.

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 × 10 −7 °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σ sensitivity for that exposure, it is cataloged as a detection.

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.
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.
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 ∼100 times.The discrepancy here is that ETAs move through the survey footprint over years rather than weeks.
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.
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 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 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 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.

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 (∼ 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.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.

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.
Figure 16 shows U (H) for the same four stability regimes shown previously.In Figure 16, there are por- 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).
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.

CONCLUSION
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.

ETA spatial distributions
When viewing L4 ETAs from Earth, the highest onsky 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 ±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.

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.

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, Ctype 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.

Figure 1 .
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°o ffset from L4.The bin amounts are scaled to the bin with the highest density.

Figure 2 .
Figure2.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 Figure1.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 Figure1.The peak of the flux-weighted density map is located 6.5°Sun-ward from L4.

Figure 3 .
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 .
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°m easured counter-clockwise from the location of Earth in the figure.Binning is the same process as done for Figure 1

Figure 5 .
Figure5.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°.

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

Figure 7 .
Figure 7.The cumulative # of unique ETA detections made by each survey for the fiducial simulation.

Figure 8 .
Figure 8.The normalized number of detections by month made by each survey.

Figure 9 .
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.

Figure 10 .
Figure 10.Asteroid recovery by each survey.

Figure 12 .
Figure 12.Total detections of ETAs by apparent magnitude.

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

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

Figure 15 .
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 Figure1.The Sun is fixed at (0°,0°), and L4 is at (0°,60°), denoted by a red dot.3) R.A. vs. Dec.

Figure 16 .
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.tions 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.

Figure 17 .
Figure 17.Cumulative upper limit on ETA populations for the three surveys.The dotted lines show the 68% CI.

Figure 18 .
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.

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