A Search for L4 Earth Trojan Asteroids Using a Novel Track-Before-Detect Multi-Epoch Pipeline
Noah Lifset, Nathan Golovich, Eric Green, Robert Armstrong, Travis Yeager
DDraft version February 19, 2021
Typeset using L A TEX twocolumn style in AASTeX63
A Search for L4 Earth Trojan Asteroids Using a Novel Track-Before-Detect Multi-Epoch Pipeline
Noah Lifset , Nathan Golovich , Eric Green , Robert Armstrong , and Travis Yeager Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, CA 94550, USA (Received 2021-2-17; Revised XXX; Accepted XXX)
Submitted to ApJABSTRACTEarth Trojan Asteroids are an important but elusive population that co-orbit with Earth at the L4and L5 Lagrange points. There is only one known, but a large population is theoretically stable andcould provide insight into our solar system’s past and present as well as planetary defense. In thispaper, we present the results of an Earth Trojan survey that uses a novel shift-and-stack detectionmethod on two nights of data from the Dark Energy Camera. We find no new Earth Trojan Asteroids.We calculate an upper limit on the population that is consistent with previous searches despite muchless sky coverage. Additionally, we elaborate on previous upper limit calculations using current asteroidpopulation statistics and an extensive asteroid simulation to provide the most up to date populationconstraints. We find an L4 Earth Trojan population of N ET < H = 14, N ET < H = 16,and N ET <
642 for H = 22. Keywords:
Asteroids, Small solar system bodies, Computational methods INTRODUCTION1.1.
Earth Trojan Asteroids
Earth Trojan Asteroids (ETAs) are an important pop-ulation of inner Solar System asteroids that could unlocka better understanding into the development and cur-rent epoch of the Solar System. They co-orbit the Sunwith the Earth at the forth and fifth Lagrange points (L4and L5) as solutions to the reduced three-body problem(see e.g. Szebehely 1967). ETAs have been estimated tobe stable on the order of 100 Myr to 1 Gyr with somespecific orbital regimes offering stability of up to the ageof the Earth ( ´Cuk et al. 2012; Marzari & Scholl 2013;Zhou et al. 2019). The oldest and most stable ETAs,referred to as primordial ETAs, could be remnants fromthe protoplanetary disk, and characterizing them couldinform our understanding of the primordial Solar Systemand its development. Studying ETAs can also provideinsight on the Earth and Moon. A leading lunar for-mation theory is the giant impact hypothesis betweenproto-Earth and Theia, which may have been a Mars-
Corresponding author: Noah [email protected] sized ETA that outgrew the stability of the L4 Lagrangepoint (Belbruno & Gott 2005). If proven, this wouldlend substantial credence to primordial ETAs. Alongwith near-Earth steroids, ETAs could explain the asym-metry in lunar impact craters between the leading andtrailing hemispheres (Ito & Malhotra 2010). ETAs arealso prime targets for space missions, because they havea low relative velocity (Malhotra 2019). This also meansthat if they are destabilized by resonances or other per-turbations, they could endanger life on Earth.While Trojan asteroids have been discovered for Venus(de la Fuente Marcos & de la Fuente Marcos 2014),Mars (de La Fuente Marcos & de La Fuente Marcos2013), Jupiter (Yoshida & Terai 2017), Uranus (de laFuente Marcos & de la Fuente Marcos 2015) and Nep-tune (Almeida et al. 2009) only one ETA has been foundto date. This is likely due to the difficulty in detectingthem as they spend much of their trajectories near theSun on the sky. ETAs co-orbit with Earth 60 degreesahead or behind (on average), corresponding to L4 andL5 respectively. As a result of this geometry, groundbased observations are limited to mostly twilight obser-vations and only short dark-sky observations dependingon the season (Whiteley & Tholen 1998). The geometryalso ensures a poor reflectance angle with the Sun, so a r X i v : . [ a s t r o - ph . E P ] F e b Lifset et al. any objects in these orbits are significantly fainter thananalogous asteroids at opposition.The only ETA found to date, 2010 TK7, was aserendipitous detection by WISE (Connors et al. 2000)at a portion of the orbit that brought it close to Earthwith a large solar elongation. Subsequent follow up sug-gests it was recently captured and is relatively unsta-ble. Its lifetime is estimated to be ∼ Previous Searches
There have been four previous searches for ETA pop-ulations. All of them found none, but three have offereda population constraint. Whiteley & Tholen (1998)used a ground based survey to search L5. They covered0.35 square degrees and set a constraint of 3 objectsper square degree at r = 22 .
8, corresponding to 350m(175m) for C (S) class asteroids. Cambioni et al. (2018)searched L4 using the OSIRIS-REX spacecraft on a flybyon its way to the asteroid Bennu. They covered nine 16deg patches and calculated an upper limit of 73 ±
23 at H = 20 .
5, which corresponds to 470/210m for C/S classasteroids. They also applied their analysis to the datafrom Whiteley & Tholen (1998), resulting in an upperlimit of 604 ±
358 at H = 21 . r = 22 .
8) and 194 ± H = 20 .
5. The Hayabusa2 spacecraft searched L5 onits way to the Ryugu asteroid (Yoshikawa et al. 2018),but published no ETA population constraint analy-sis beyond a null detection. Markwardt et al. (2020)searched for ETAs at L5 using the Dark Energy Cam-era. They covered 24 square degrees and calculated anupper limit of N ET < H = 15 . N ET < − − H = 19 .
7, and N ET <
97 at H = 20 .
4; this isthe strictest constraint on the ETA population to date.In this paper we present the results of a new searchfor ETAs at the Earth-Sun L4 point and the new pop-ulation upper-limit calculated from a null result. In § §
3, wedescribe the detection pipeline and the ETA simulationwe used to determine our upper limits. In §
4, we presentthe results of our detection pipeline and the upper limitcalculation. In §
5, we compare our results to previous searches and provide a meta-analysis of ETA populationupper limits across the literature. DATA2.1.
DECam Observations
Our survey consisted of roughly 90 minutes of dataon both 2019 Nov 20 and 21 on the Blanco telescope us-ing the Dark Energy Camera (DECam Honscheid & De-Poy 2008). DECam has 60.5 active 4k ×
2k CCDs with0 (cid:48)(cid:48) . , pixel − . We observed in the r-band, choosingslightly lower solar reflectance from potential ETAs forlower sky background during twilight. We unfortunatelyfailed to include dithers in the observing scripts. How-ever, since our detection pipeline stacks signal across allexposures, we are still sensitive to ETAs that traversechip-gaps during our survey time.Night 1 and 2 consist of 78 and 77 exposures, respec-tively, with all exposures from a given night stacked onthe same position. The exposures were all 40 seconds,with roughly 30 seconds between exposures for readout.We covered an area of 5.72 deg located near the L4Lagrange point (see Figure 1). This effective area is es-timated by a computing the area of a circle of radius1.1 ◦ , while accounting for the overlap in pointing be-tween the two nights. We selected a location slightlyeast of the Lagrange point because the peak of the onsky spatial distribution is located there as a result ofthe angled observation geometry and the asymmetry oftadpole Trojan orbits. ETAs are expected to have amotion of ∼ ◦ day − . This allows for consecutive-nightdetections for typical ETA trajectories on sky.2.2. Reduction Pipeline
We obtained our exposures after standard process-ing by the DECam Community Pipeline (Valdes et al.2014). We then used the LSST software package (Boschet al. 2019) to produce difference images. This can besplit into three different pieces: making calibrated sin-gle epoch catalogs and images, creating a static tem-plate image, and subtracting the template from the sin-gle epoch images to produce difference images.The LSST stack single epoch processing includes: themasking of cosmic rays, measuring the PSF, detectingobjects, deblending and measuring individual sourcesand calibrating the astrometry and photometry. A de-tailed explanation of each of these steps is explained inthe HSC data release papers (Aihara et al. 2018, 2019).For photometric calibration, we used the Pan-STARRScatalog (Flewelling et al. 2016). https://pipelines.lsst.io rack Before Detect Earth Trojan Survey Figure 1.
The data set depth (number of frames) in R.A.and Dec. We cut the entire region into “patches” with theLSST data management pipeline, followed by creation of amaster pixel grid. All analyses occurred on this grid.
To create the template images we first interpolatedthe data onto a tangent plane projection centered onthe average pointing. The total observed area was di-vided into a regular grid with each grid 4000 × ∼ (cid:48) on a side (see Figure 1). For the tem-plates we selected all images with PSF full-width at halfmaximum seeing < . (cid:48)(cid:48) . For each grid the static skyis constructed following the procedure outlined in (Ai-hara et al. 2018, 2019). Briefly, the procedure builds atwo-sigma clipped coadd to construct a static image ofthe sky. We subtract this coadd from each individualimage and identify variable sources. Those detectionsthat are truly variable will only appear in a small sub-set of the single epoch visits. The variable sources arethen masked and a coadd image is created by taking themean of all the images.We use the Alard-Lupton algorithm ( ? ) as imple-mented in the LSST stack to create difference images.This procedure estimates a convolution kernel which,when convolved with the template, matches the PSF ofthe template with that of the science image by minimiz-ing the mean squared difference between the templateand science image. The Alard-Lupton procedure useslinear basis functions, with potentially spatially-varyinglinear coefficients, to model the matching kernel whichcan flexibly account for spatially-varying differences inPSFs between the two images. The algorithm has theadvantage that it does not require direct measurementof the images’ PSFs. Instead it only needs to modelthe differential matching kernel in order to obtain anoptimal subtraction. After examining the difference images, we found thatthere remained a large number of artifacts due to prob-lems in the difference imaging algorithm. A majority ofthese objects were bright stars, but even fainter starscaused problems for our analysis. We tried to removethese issues by masking out known objects. We selectedobjects from PAN-STARS with m r <
20 and maskedthem using the procedure described in Coupon et al.(2018). For each star the mask was composed of amagnitude-dependent circle for the star and a rectan-gle for the bleed trail. We visually tuned the size of thecircle and rectangle to match our data. We split thedata into two class above and below 14th magnitude.The size of the radii in arcseconds are:radius = e − m r / . , for 14 < m r < e − m r / . , for m r < . (1)The length and width of the rectangles are 1.5 × < m r <
20 and6 × m r <
14. The longer rect-angles were necessary to remove some long bleed trailsfor bright stars. METHODSWe employ a shift-and-stack method to detect aster-oids in the data. By combining the signal from multi-ple images, we can reduce the noise floor and increasesensitivity. We use a novel shift-and-stack method de-veloped by (Golovich et al. 2020) that can be describedas track-before-detect. This works by generating a largenumber of random potential asteroid trajectories withinthe search priors, calculating the signal-to-noise-ratio(SNR) for each trajectory in the data, and saving theinformation of potential detections if the SNR meetsa given threshold. This process is then repeated un-til the priors have been completely filled through therandom trajectory generation. Traditional track-before-detect methods will usually stack entire images on topof each other according to a selected angle or proper mo-tion (e.g. Heinze et al. 2015). This can work just as welland will save some computational time since all paral-lel trajectories are simultaneously sampled by detectingon the shifted stack. However, this is restricted to lin-ear source motion, because different proper motions cannot be applied to different parts of the same stackedimages. Track-before-detect allows for non-linear tra-jectories at the (substantial) cost of sampling trajecto-ries individually, because all that is needed is the po-sitions of intersections between a given trajectory andimages in the data set. In this paper, we use linear tra-jectories, which were sufficient to generate competitive
Lifset et al.
Figure 2.
A flow chart demonstrating the detection pipelinefrom start to finish. Ovals represent filters and boxes repre-sent calculations or trajectory interactions. Numbers of po-tential detections after each step are shown on the right forboth nights. The final potential detection lists are manuallysorted to produce final detection lists. Differences betweenthe two nights are due to variation of difference imaging ar-tifacts, sky noise variation, and inherent randomness. results in an ETA constraint; however, we implementour track-before-detect pipeline as a second step (follow-ing Golovich et al. 2020) toward non-linear trajectoriesand applications to larger data sets spanning much moretime and sky coverage.3.1.
Detection Pipeline
The detection pipeline is shown as a flow chart in Fig-ure 2 with the number of potential detections shown ateach step for both nights. We run our detection algo-rithm separately on each night with the same priors, be-cause ETAs can not be expected to move linearly on timescales that would include both nights of data. We rep-resent each night’s survey data in memory as a virtual,linear, three-dimensional coordinate space that encom-passes the minimum bounding box to contain the sur-vey’s image data. In this space, the X- and Y- axes cor-respond to individual pixels in the image data while theZ-axis represents time in seconds. This also allows us tovirtually stitch together the images from each exposureby placing the individual image data at the correct posi-tions within this three-dimensional data space. Each ex- posure represents 40 seconds of real time, so pixels maybe imagined as rectangular prisms with cross-sectionalarea equal to the pixel area, and height proportional tothe exposure timeThe pipeline begins with trajectory generation withinthe defined priors. Each trajectory is defined by 5numbers: X-intercept, Y-intercept, Z-intercept/time, X-slope, and Y-slope. The X- and Y-intercepts are thepixel position of the trajectory at that Z-intercept/timeand the X- and Y-slopes are the corresponding velocitiesin pixels s − . In trajectory generation, the initial valuesof the X-, Y- and Z-intercepts were chosen by randomlyselecting a non-NaN pixel in the entire night’s data andpropagating back to the beginning of the first exposureto calculate intercept values. The slopes were randomlygenerated from a uniform rectangle in ecliptic coordi-nates: -40 to 40 (cid:48)(cid:48) hour − in longitude and 110 to 190 (cid:48)(cid:48) hour − in latitude. This range was chosen based onWiegert et al. (2000) to cover nearly all potentially ob-servable ETAs and aligns with the range used by Mark-wardt et al. (2020) (0.75–1.25 ◦ day − ). We generated8 × trajectories for each night in order to fill theparameterized space. This number was chosen empiri-cally as the recovery rate of injected asteroids no longerimproved with additional trials.3.1.1. Pixel Values and Memory Management
Once the trajectories are generated, we calculate theirintersections with each of the images in the data set andpull the pixel values for each of those points. However,the large number of trajectories and the full volumeof image data make this process complicated in prac-tice. We first split the total number of trajectories intoequally-sized work units to be distributed among workernodes on our computing cluster. These work units arefurther divided into a series of batches by the workers.For each night of data, we tested approximately 824 bil-lion trajectories across 768 work units. Each work unitevaluated approximately 536 million trajectories thatwas further subdivided into 8 batches of approximately67 million individual trajectories drawn from the motionmodel priors.Each batch is processed by sequentially streamingthrough all of the exposures in the survey data and col-lecting the information needed from each image. Forevery intersection point in an exposure, we calculate aweighted sum of the image pixels and store the localflux, the number of non-NaN pixels intersected, and lo-cal SNR (see § rack Before Detect Earth Trojan Survey .3.1.2. SNR Calculation
Once the requisite pixels are pulled from the data,the SNR is calculated with a weighted sum on the opti-mal signal matched filter using a two-dimensional Gaus-sian profile convolved with a line segment equal to thelength of the expected streak, which we can computefrom the sampled proper motion for each trial trajec-tory and the exposure time. The flux and noise in eachpixel were summed across all images that a given sourceintersected. Within one image intersection, we measurethe flux F by calculating a weighted sum over the pixels( p ) inside the streak: F = (cid:88) i w i p i . (2)The weights ( w ) are determined from the profile, whichis the convolution of a line segment with a two-dimensional Gaussian. The weights are normalized suchthat their sum equals the effective number of pixels inthe streak: n eff = (cid:88) i w i . (3)The effective number of pixels is calculated by integrat-ing the flux profile: n − eff = (cid:90) dxdy P ( x, y ) . (4) https://computing.llnl.gov/computers/catalyst Specifically, we consider a streak with length L orientedalong the x-axis, width σ y , flux per unit length l , andcenter (0,0), convolved with a symmetric bivariate Gaus-sian PSF of width σ π : P ( x, y ) = l L (cid:113) π ( σ π + σ y ) exp (cid:18) − y / σ π + σ y (cid:19) × (cid:32) erf (cid:32) L − x (cid:112) σ π (cid:33) + erf (cid:32) L + 2 x (cid:112) σ π (cid:33)(cid:33) . (5)Assuming that the object is unresolved, σ π >> σ y , thissimplifies accordingly: P ( x, y ) = l L √ πσ exp (cid:18) − y σ (cid:19) × (cid:18) erf (cid:18) L − x √ σ (cid:19) + erf (cid:18) L + 2 x √ σ (cid:19)(cid:19) , (6)where we have defined σ ≡ σ π .We choose an area surrounding the streak to sum overthat balances capturing as much signal as possible withlimiting computation time. We integrate flux to 3 σ as-suming a Gaussian point spread function away from thestreak’s ridge, which encompasses over 99% of the sig-nal.Sources of noise include sky background Poisson noise,read noise, dark current, and Poisson noise from thesignal itself. The total noise N includes each of thesebut is dominated by the sources of Poisson noise, whichadd in quadrature. Thus we have: SN R = F √ F + N . (7)The background noise is calculated as the standard de-viation of pixel values in the difference image. The ex-pected value per pixel will then be the square of thisnumber, and thus we get the total background signal as N = σ n eff . (8)From Equations 2-6, the flux from the asteroid is givenby F = (cid:80) i ( l i p i ) (cid:80) i l i (cid:80) i l i , (9)where the sum over pixels approximates the integral.There is a chance that a valid detection might overlapon one image with an image artifact like a cosmic ray.To compensate for this, we store an initial light curve foreach trajectory consisting of the local SNR value on eachintersected image. We then apply an outlier rejectionmethod, consisting of a rolling median sigma clip, to thistemporary light curve and ignore those clipped images Lifset et al. when calculating the final SNR for the trajectory. Westore information on any trajectory that has an SNRgreater than 5 and intersects at least 5 images. Mostof these (see Figure 2) are false positives or duplicatedetections of the same object that will later be filteredout. 3.1.3.
Detection Post-processing
Once we have initial detection lists for each night, weapply a filter to remove duplicate detections of the sameobjects. This is common for bright asteroids, because atrajectory that is only partially aligned with the aster-oid will still surpass the SNR threshold. We use a k-dtree to find nearest neighbors for all the trajectories inthe four dimensional space of ( x, y ) position and propermotion, where we convert the proper motions to spatialcoordinates by multiplying by the total integration timeof the data for each night. Using this four-distance met-ric, we consider all trajectories within (cid:15) = 5 (cid:48)(cid:48) of othersduplicates. To be explicit, we require: (cid:113) ( | (cid:126)x ,i − (cid:126)x ,j | ) + (∆ t | (cid:126)π i − (cid:126)π j | ) < (cid:15), (10)where (cid:126)x ,ij are the position of two trajectories at thefirst exposure, ∆ t is the time between the last and firstexposure, and (cid:126)π ij are the proper motions. We keep thetrajectory with the highest SNR. We chose (cid:15) = 5 (cid:48)(cid:48) asa balance between removing a meaningful number ofduplicates and being small enough to have almost nochance of removing any two different objects that areclose by coincidence. This cut is small enough that manyduplicate trajectories will still pass through the filter,which will be addressed by a stronger duplicate filter atthe end of the pipeline.We then refine each trajectory with a Scipy (Virtanenet al. 2020) optimizer using the Nelder-Mead methodto maximize global SNR. We save small image cutoutsof each trajectory-image intersection and use the sameSNR calculation as above on these image cutouts as anoptimization function. We then take these refined tra-jectories and apply a few stronger filters to remove po-tential false positives.First, we apply a series of filters on the median stackimage. This median image is calculated by taking themedian across the image cutouts for each trajectory. Inthis median image, a trajectory that perfectly alignswith a source would appear as a centered single PSF(or short streak for fast moving sources). We calculatethe SNR of this median stack using a the same weightedsum as before and we require
SN R med >
7. Next wecalculate the the χ of the center of the median imageand the χ of an annulus region around it. To be ex- plicit, we calculate: χ = 1 N pix (cid:88) p i N i (11)where the center consists of pixels within 3 σ of the cen-ter and the annulus consists of pixels between 5 σ − σ ,where σ is with width of the Gaussian PSF. We requirethat χ annulus < χ center < χ annulus −
1. The ideabehind this filter is to discard false positives based onshape. False positives that result from difference im-age artifacts usually have signal in the annulus regionas well as the center. The exact cut was determined em-pirically by calculating the χ values for a set of knowninjected asteroids and setting limits that would balancediscarding the false positives with not discarding anyvalid detections.Next, we apply a “near-hit” filter that attempts to cutout any trajectory that only partially overlaps with areal object in the data. These are often duplicate detec-tions of bright objects as previously discussed, but aretoo far in 4-distance from the center to be filtered outthat way. The near-hit filter is carried out by analyzingthe distribution of local SNR values in the light curve.For a good hit, most of the local SNR values would dis-tribute around a mean value due to Poisson noise in theflux and the sky background as well as intrinsic variabil-ity due to rotation of the asteroid. For a near-hit, theportion of the trajectory that misses the true asteroidmotion would exhibit near-zero SNR. To differentiatebetween these two, we rank-order the local SNR valuesinto an array and compared the standard deviation ofthe highest values of the distribution (fourth quartile) tothe standard deviation of the lowest (first through thirdquartile). We required that the fourth quartile variancebe less than the first through third quartile variance,thus requiring that most of the high value local SNRsare clumped around a value indicative of a detectionon all of the frames that intersect the trajectory. Wedemonstrate asteroids that both pass and fail this cutin Figure 3. In the case where all values are clumpedaround zero, and this filter is passed, the overall SNRwould likely be below our threshold in the first place.We also apply a similar filter where we take each trajec-tory and randomly split its light curve in two. We thenrequire that each half individually must reach the SNRcut of √ to check that flux is roughly split between tworandom halves of the light curve.After false positive filtering, we apply another strongerduplicate filter with a 4-distance cutoff of 30 (cid:48)(cid:48) . Thiswas chosen because it is the same 4-distance used formatching detections with known objects like MPC cat-alog asteroids or injected false asteroids (see § rack Before Detect Earth Trojan Survey Figure 3.
Example distributions of Local SNR values offrames hit for detection trajectories that are bad (left panel)and good (right panel) approximations of the true trajectory.The near hit filter compares the standard deviation of thefirst through third quartile (blue bins) with that of the fourthquartile (orange bins). The near-hit mostly has local SNRsaround 0 and some reaching higher corresponding to whenthe trail trajectory overlaps the flux from the asteroid, andso fails the test. The right panel shows a good hit, whichhas mostly non-zero local SNRs distributed around a typicalvalue related to the brightness of the asteroid. final duplicate filter reduces the potential detection listto a small number that can be manually sorted and ul-timately made into a final detection list.3.2.
False Source Injection
In order to test completeness, we create an identicaldata set with false asteroids injected during data reduc-tion. Roughly 4000 false asteroids were injected sepa-rately in each night. They were injected as linear tra-jectories with slopes randomly selected from the samepriors used in detection. Their X-, Y-, and Z-interceptswere generated by randomly selecting a valid right as-cension, declination, and time within the bounds of thenight’s data. They had apparent magnitudes randomlyselected between 17 and 26. Any potential injectionsthat intersect fewer than five images were not consid-ered in the completeness calculation.We determined which detections corresponded to in-jected asteroids using the 4-distance metric defined inEquation 10. A given injected asteroid was consideredmatched with a detection if there was one within a 4-distance of 30 (cid:48)(cid:48) and the closest one was selected if therewere multiple. We chose this cutoff empirically. Wefound that 30 (cid:48)(cid:48) was large enough to include nearly allvalid detections including those that did not align per-fectly. It was small enough to rarely include two actuallyseparate injected asteroids that were by chance overlap-ping (see Figure 4). Once matches are made we are
Figure 4.
Distribution of 4-distances from each injected as-teroid to the closest other injected asteroid. We show thisfor each asteroid instead of each closest pair, so all numbersare double counted. A cutoff of 30 (cid:48)(cid:48) was used despite the po-tential to rule out two genuinely separate injected asteroids,because it was needed to include nearly aligned detectionsin order to reduce overall potential detections more signifi-cantly. able to measure completeness as a function of apparentmagnitude. 3.3.
ETA Simulation
In order to calculate a constraint on the ETA popu-lation, we determine the fraction of a theoretical ETApopulation that our observations cover. This is done bysimulating the orbits of millions of asteroids over 10 years to determine the spatial distribution of long-termstable ETAs. One of the main improvements this papermakes over previous ETA population calculations is theuse a large scale simulation in determining the distribu-tion of stable ETAs for the constraint calculation.The simulation was run with the Rebound
N-body in-tegrator software using a IAS15 integrator with adaptivetime steps (Rein & Liu 2012; Rein & Spiegel 2015). Weinitialized roughly 4 million point particles in a simula-tion including the Sun, the Earth, and Jupiter. Incli-nations were randomly selected from a Gaussian dis-tribution centered on 0 with a standard deviation of10 degrees. Eccentricities were randomly selected froma Gaussian distribution centered on 0 with a standarddeviation of 0.1 (absolute value taken after sampling).Semi-major axes were randomly selected from a Gaus-sian distribution centered on 1 AU with a standard de-viation of 0.1 AU. All other orbital elements were ran-domly selected from a uniform distribution from 0 to2 π . If at any point in the simulation, a point particlemoved outside the half-sphere on the L4 side of the or-bit relative to the plane perpendicular to the eclipticthat contains the Earth–L3 line, it was no longer con- Lifset et al.
Figure 5.
Spatial distribution of stable ETA population.Points are sampled every 10 years for 1 million years. Theblack dot represents the Lagrange point and red circles thetwo nights of observations. Positions at each time step areincluded in order to show the time average position of ETAs. sidered a stable ETA and was discarded. This criteriawas focused on finding long term stable ETAs includ-ing large tadpole orbits but not horseshoe orbits thatcontinue beyond L3 to the other half-sphere. Of the ini-tial 4 million point particles with Earth-like orbits, thesimulation ended with 145,000 stable ETAs under theinfluence of the Earth and the Sun. Positions are sam-pled every 10 years and included together to show timeaverage positions. Their geocentric ecliptic coordinatesare calculated and shown in Figure 5. RESULTS4.1.
Detections
We found no ETA candidates in either night of obser-vations. Our priors were strict and focused on potentialETAs, so we expected to detect few objects in general.We detected 55 objects: 18 on night 1 and 37 on night 2.All except one were known objects in the MPC catalog.The single new detection was found to have a propermotion of 0.58 ◦ day − , which was outside the expectedrange for ETAs. It is likely that we detected this object,as well as the others, on the edge of our priors and therefinement process aligned it more properly outside thepriors. Its median image and light curve are shown inFigure 6. We detected four Mars-orbit crossing Amortype asteroids and one Earth-orbit crossing Apollo typeasteroid. Our ability to detect an Apollo asteroid andfour Amor asteroids is a testament to our pipeline’s abil-ity to detect asteroids of interest in the same areas thatan ETA would exist. The pipeline can also be utilized todetect these orbital regimes preferentially with a differ-ent choice in priors. All detected asteroid information isshown in the appendix in Table 2 and Table 3. Figure 6.
Left: Median stack image of the new detectedasteroid. The increased noise at the top is because the CCDedge is nearly aligned with the top of this cutout. Right:light curve of new detected asteroid with local SNR as afunction of image intersected. The drop at the end is from astar mask that overlaps with the trajectory.
ETA Upper Limit
We convert from PAN-STARR’s calibrated r-bandmagnitude to apparent magnitude V: V = m r + 0 . , (12)which assumes ETAs are most similar to S-class aster-oids with albedo 0.2 (Chesley & Veres 2017). We thenconvert to absolute magnitude H: H = V + 2 . ( φ ) , (13)where the phase integral φ is given by φ ( G = 0 .
15) = 0 .
85 exp (cid:16) − .
332 tan ( α/ . (cid:17) +0 .
15 exp (cid:16) − .
862 tan ( α/ . (cid:17) , (14)and α is the Sun-asteroid-Earth angle (Dymock 2007).Here we assume a distance of 1 AU and angle of 60 ◦ based on the L4 Lagrange point. Our recovery rate is >
50% for
V < .
92 and
H < .
77. Our recoveryrate as a function of absolute magnitude H is shown inFigure 7.We calculate our simulated ETA population coverageusing the time average position of stable ETAs over thecourse of the simulation shown in Figure 5. We findthat our observations cover 0 .
91% of the ETA popula-tion. There is some sensitivity to simulation initial con-ditions, but we find this distribution to be most repre-sentative of a potential ETA population given the largenumber of asteroids we injected into our simulation andthe long timescale of our simulation during which theETAs come into equilibrium with the Earth and Sun.Our coverage is substantially lower than previous ETAsearches, because we assume a much broader stable ETApopulation based on simulation results. Even with vary-ing simulation initial conditions, our population is stillsubstantially broader. We discuss coverage comparisons rack Before Detect Earth Trojan Survey Figure 7.
Recovery rate of injected asteroids as a func-tion of absolute magnitude H. The red dotted line shows50% recovery and thus overall sensitivity, which was 21 . and the significance of our simulation in the constraintcalculation in more detail in § σ , by Poisson statistics. Our upper limit is then U ( H ) = 3 R ( H ) ∗ C , (15)where R ( H ) is our recovery as a function of absolutemagnitude H and C is our coverage percentage. Ourcalculated upper limit is shown in Figure 8. This cal-culation is for the L4 Lagrange point only. We assumethat L5 has the same upper limit by symmetry and so wealso freely compare with previous searches of L5. Theupper limit curve can be broken into two areas. The flatupper limit on the left and middle is a result of our max-imum recovery rate through an H magnitude of 20 . U ( H ).If we assume that ETAs follow a similar power lawH-distribution as Near Earth Objects (NEOs), we cancalculate a more stringent and accurate U ( H ). We makeuse of the H distribution for NEOs: N ( < H ) = A ∗ αH , (16) Figure 8.
Upper limit of the ETA population at L4 directlycalculated from our recovery rate. The flat area on the leftis a result of maximum recovery with small coverage. Thesloped area on the right is where our recovery starts to drop. where α = 0 . ± .
02 for 13 < H < α = 0 . ± .
01 for 16 < H <
22, and α = 0 . ± .
03 for
H >
22 (Schunov´a-Lilly et al. 2017). This is consistent withHarris & D’Abramo (2015); Tricarico (2017); Morbidelliet al. (2020). We take our strongest upper limit resultof N ( < H = 21 . .
17 and extrapolate usingEquation 4.2 (see Figure 12).
Table 1 . Constraint Analysis ResultsStudy Narrow ETA Population Narrow ETA Population Broad ETA Population Broad ETA PopulationMarkwardt H Distribution Current H Distribution Markwardt H Distribution Current H DistributionThis Study 8.74 22.40 54.81 140.45Markwardt 23.60 21.77 134.02 123.62Whiteley and Tholen 108.30 242.68 1174.94 2632.78
Table 1 continued Lifset et al.
Table 1 (continued)
Study Narrow ETA Population Narrow ETA Population Broad ETA Population Broad ETA PopulationMarkwardt H Distribution Current H Distribution Markwardt H Distribution Current H DistributionCambioni 40.45 49.93Comparisons of ETA population constraints at a single Lagrange point at H = 20 (297 meters at 0.2 albedo). We comparethe constraints set by our analysis, Markwardt et al. (2020), Whiteley & Tholen (1998), and Cambioni et al. (2018).5. DISCUSSION5.1.
Comparison to Previous Searches
The two major differences between our upper limitcalculation and previous searches are different H distri-butions and a broader ETA population used for deter-mining coverage. The H distribution used in Cambioniet al. (2018); Markwardt et al. (2020) is a power law with α = 0 . α = 0 . − .
7, respectively. Our H distri-bution is a power law with α = 0 .
48 for 13 < H < α = 0 .
33 for 16 < H <
22, and α = 0 .
62 for
H > C = 5 .
72 for the narrowETA population. We estimate the coverage of Whiteley& Tholen (1998); Markwardt et al. (2020) of the broadETA population using their given observational informa-tion. Unfortunately, Cambioni et al. (2018) provides toolimited observational information, but their constraint isweaker than ours and that of Markwardt et al. (2020).We repeat our analysis for four situations, which areshown in Figures 9, 10, 11, and 12 with results for H = 20 shown in Table 1. Here we assume the ETA pop-ulation at L4 and L5 are identical. We consider 12 themost accurate and up to date constraint on the EarthTrojan Asteroid population. When using our more ac-curate H distribution, our upper limit is consistent withMarkwardt et al. (2020). When using the H distribu-tion of Markwardt et al. (2020), we have a stronger re-sult. The broad ETA population significantly increasesthe upper limit calculated by all searches. The upperlimit values depend heavily on assumed ETA popula-tion distribution, and could vary if a different popula-tion distribution is used. We are currently preparing afull analysis of different ETA simulations and potentialpopulation distributions (Yeager et al. in prep).5.2. Outlook for ETA Surveys
In this paper we demonstrate that stacking signalacross multiple epochs is an effective way to increase sur-vey sensitivity for ETAs. This is not surprising, as thismethod has been used for many populations of aster-oids. However, this technique has important limitationswhen applied to ETAs. All inner solar system asteroidsenter the regime of non-linear motion within approxi-mately a single night (see Figure 4 of Heinze et al. 2015,for a clever demonstration of this). This limits the timeover which signal may be stacked to increase sensitivity.We have employed this fact to reach comparable depthsof Markwardt et al. (2020) despite a quarter the on skycoverage.Where our method could excel beyond this limitationis in its track-before-detect framework, which enablesdetection even on non-linear trajectories. This requires rack Before Detect Earth Trojan Survey Figure 9.
Upper limit of the ETA population at L4 us-ing a narrow ETA population and the H-mag distribution ofMarkwardt et al. (2020). Results are shown for this paperand the other three from the literature (Whiteley & Tholen1998; Cambioni et al. 2018; Markwardt et al. 2020).
Figure 10.
Upper limit of the ETA population at L4 usinga narrow ETA population and our up-to-date H-distribution.Results are shown for this paper and the other three from theliterature (Whiteley & Tholen 1998; Cambioni et al. 2018;Markwardt et al. 2020). the parameterization of the motion model in order totrack for arbitrarily long timescales and requires theability to quickly generate large numbers of trajectories.However, ETAs move under comparable influence of theSun and the Earth, and thus their trajectories are notparameterizable by Keplerian orbital elements. Thereis promise in several machine learning methods (con-volutional neural networks, generative adversarial net-work, and spatial statistical methods involving Gaussianprocesses) to develop surrogate models to quickly gen-erate ETA trajectories to feed our track-before-detectpipeline. However, this is beyond the scope of this pa-
Figure 11.
Upper limit of the ETA population at L4 us-ing our broad ETA population and the H-mag distributionfrom Markwardt et al. (2020). Results are shown for this pa-per along with Whiteley & Tholen (1998); Markwardt et al.(2020).
Figure 12.
Upper limit of the ETA population at L4using our broad ETA population and our up-to-date H-distribution. Results are shown for this paper, along withWhiteley & Tholen (1998); Markwardt et al. (2020). Weconsider this to be the most current and accurate represen-tation of the ETA population. We find an L4 Earth TrojanPopulation of N ET < H = 14 . N ET < .
72 for H = 16, and N ET <
642 for H = 22. per. We are eager to learn of new developments on thisavenue as large scale twilight surveys could unlock in-credible potential to detect ETAs. Until then, shift-and-stack methods on single night data sets carried out overincreasingly more nights with larger telescopes and cam-eras with larger fields of view and shorter overheads willbe the best we can do.2 Lifset et al.
ACKNOWLEDGMENTSThis work was performed under the auspices of theU.S. Department of Energy by Lawrence Livermore Na-tional Laboratory under Contract DE-AC52-07NA27344and was supported by the LLNL-LDRD Program underProjects 17-ERD-120, 19-SI-004, and 20-ER-025. Thiswork was based on observations obtained at Cerro TololoInter-American Observatory a division of the NationalOptical Astronomy Observatories, which is operated bythe Association of Universities for Research in Astron-omy, Inc. under cooperative agreement with the Na-tional Science Foundation.
Facilities:
CTIO(Blanco/DECam)APPENDIX
Table 2 . Night One DetectionsMPC Designation V R.A. Dec. π α cos δ a π δ Observation Time b Magnitude ◦ ◦ ◦ day − ◦ day − MJD(1) (2) (3) (4) (5) (6) (7)28445 19.2 299.00217 -20.80504 0.4 0.046 58807.0124656588 18.2 299.6919 -20.43834 0.35 0.059 58807.01246525990 19.0 300.15762 -20.05189 0.569 0.082 58807.01326418284 18.5 300.0047 -20.69366 0.629 0.166 58807.012465382363 22.0 300.01745 -20.2974 0.626 0.119 58807.014861218860 21.6 299.60271 -19.96368 0.53 0.056 58807.012465474354 22.2 298.72499 -21.01137 0.565 0.143 58807.0516271096 18.5 299.2128 -20.24981 0.367 0.111 58807.012465359369 22.4 300.60024 -20.9554 0.643 0.106 58807.026817428658 21.1 299.39684 -20.28277 0.688 0.088 58807.012465321997 21.3 299.72667 -20.82516 0.634 0.089 58807.02203766597 20.2 299.03489 -20.69781 0.561 0.091 58807.01246530031 19.5 299.52837 -21.07535 0.493 0.139 58807.0124652763 16.6 299.08874 -19.91832 0.524 0.13 58807.022037474354 22.2 298.72479 -21.01091 0.651 0.139 58807.054028158 15.0 299.31813 -20.15988 0.313 0.061 58807.012465933 18.0 300.0497 -20.13273 0.319 0.068 58807.04206200055 20.2 300.56029 -20.24477 0.547 0.151 58807.012465 aα and δ refer to right ascension and declination, respectively. b Observation time refers to the time at which the R.A. and Dec. values were observed. These were the firstexposure the detected trajectory intersected. rack Before Detect Earth Trojan Survey Table 3 . Night Two DetectionsMPC Designation V R.A. Dec. π α cos δ π δ Observation TimeMagnitude ◦ ◦ ◦ /day ◦ /day MJD(1) (2) (3) (4) (5) (6) (7)19.7 300.81914 -20.28069 0.559 0.171 58808.0121766429 18.3 299.92129 -20.854 0.579 0.097 58808.01217619890 20.5 300.67434 -20.1136 0.41 0.184 58808.01217679637 20.2 300.10334 -19.72785 0.499 0.076 58808.014583933 18.0 300.39079 -20.0857 0.329 0.053 58808.0121766186 18.6 301.1003 -19.93354 0.404 0.102 58808.0536117316 19.9 301.4632 -21.00672 0.305 0.119 58808.01217630031 19.5 300.06307 -20.93398 0.508 0.153 58808.012176137296 20.4 301.31721 -20.52184 0.56 0.089 58808.01217633639 18.7 300.41655 -20.58602 0.346 0.014 58808.01217611464 20.9 299.81264 -20.09294 0.241 0.056 58808.012176428658 21.1 300.13527 -20.19834 0.687 0.089 58808.01217634004 18.5 301.62358 -20.22192 0.532 0.123 58808.01217633056 20.3 299.76062 -20.77199 0.359 0.113 58808.06475732002 19.7 300.33409 -20.1858 0.517 0.08 58808.012975153340 20.0 301.44698 -20.4279 0.559 0.154 58808.0321764795 18.3 300.93586 -21.21965 0.514 0.08 58808.012176183248 21.3 300.20088 -20.77253 0.493 0.069 58808.01217638240 20.2 300.57507 -20.75502 0.5 0.08 58808.01217667552 20.7 300.92415 -19.76929 0.597 0.091 58808.03692127162 18.7 300.09441 -20.9108 0.389 0.112 58808.013785245878 21.1 300.24139 -19.76707 0.485 0.086 58808.03692144031 20.8 300.9262 -19.56985 0.376 0.034 58808.01943367729 20.6 301.67756 -20.77597 0.627 0.1 58808.01297522520 19.4 300.55075 -20.0181 0.558 0.098 58808.01217640655 19.8 300.06836 -20.56463 0.508 0.081 58808.0121766588 18.2 300.06775 -20.37796 0.342 0.08 58808.01217618284 18.5 300.68213 -20.5291 0.63 0.165 58808.01217698050 20.3 301.10895 -20.68522 0.548 0.132 58808.0121766421 18.4 301.10738 -20.32112 0.38 0.065 58808.012176138278 20.3 301.56737 -20.94839 0.475 0.111 58808.012176200055 20.2 301.14832 -20.0934 0.556 0.145 58808.01217646863 20.1 301.28386 -20.2624 0.349 0.073 58808.05282437461 20.8 300.73613 -19.76455 0.556 0.099 58808.024236124069 21.1 301.09555 -21.41326 0.246 0.13 58808.05601962052 20.4 300.06279 -20.11545 0.37 0.088 58808.012176322045 21.7 300.09013 -21.20186 0.583 0.129 58808.05838 Lifset et al.