Synthetic Tracking Using ZTF Deep Drilling Data Sets

The Zwicky Transit Facility (ZTF) is a powerful time domain survey telescope with a large field of view of 47 deg 2 . We apply the synthetic tracking technique to integrate a ZTF’s deep drilling data set, which consists of 133 nominal 30 s exposure frames spanning about 1.5 hr, to search for slowly moving asteroids down to approximately 23rd magnitude. We found 1168 objects from searching 40 of the 64 CCD-quadrant subfields, each of which covers a field size of about 0.73 deg 2 . While most of the objects are in the core region of the asteroid belt, there are asteroids belonging to families of Trojan, Hilda, Hungaria, Phocaea, and near-Earth-asteroids. Such an approach is effective and productive for discovering new asteroids. Here we report the data processing and results as well as discuss a potential deep drilling operation mode using this approach for survey facilities.


Introduction
Aligning and stacking up (shift-and-add) images to improve sensitivities has been performed for some time to search for faint outer solar system objects (Tyson et al. 1992;Cochran et al. 1995;Bernstein et al. 2004). Recently, we developed the synthetic tracking (ST) technique based on the same idea to search for near-Earth-objects (NEOs) by (1) acquiring frames fast enough to avoid streaks in a single frame and (2) searching systematically for moving signals in post-processing Zhai et al. 2014) simulating tracking at velocities in a range. Sophisticated data preprocessing and massively parallel computation is the key to success. The signals are not required to be detectable in individual frame as is typically the case for detecting faint objects. Massively parallel computation is needed to do a systematic search over a velocity grid covering a velocity range of interest. Modern GPUs are typically used for the purpose Whidden et al. 2019). We adopt GPU-aided computation to carry out this search by brute force and are able to keep up the real-time processing of data from a survey. This allows us to integrate multiple frames for a long time to gain sensitivity for detecting faint asteroids without suffering trailing loss. For detecting moving objects, with ST, a small telescope can achieve sensitivity of a large telescope with sufficiently long integration time.
Synthetic tracking can be applied to different timescales depending on the rate of motion of objects of interest (Heinze et al. 2015). The rule of thumb is to set exposure time such that the fastest motion does not streak, i.e., the motion is no more than the size of the point-spread-function (PSF) measured as full width at half maximum (FWHM) or one pixel depending on which one is larger. For searching faint satellites and NEOs, modern CMOS cameras with low read noise (typically 1-2e − ) are required to take frames fast enough to "freeze" these objects in a single exposure, yet still limited by the sky background.
The Zwicky Transit Facility (ZTF) performs time-domain survey with a 47 square degree field of view covered by a mosaic of 16 CCDs mounted on the Samuel Oschin 48 inch Schmidt telescope Graham et al. 2019). One of the goals of ZTF is to survey asteroids, which is important for planet protection, solar system formation study, supplying mission targets, as well as asteroid-mining. The nominal ZTF data processing finds asteroids as significant signals in a single frame with an exposure time of 30 s taken by CCDs (Masci et al. 2018;Mahabal et al. 2019). Here we present results from applying ST to a ZTF deep drilling data set consisting of 133 30 s exposure frames over about 1.5 hr observing the same field. Applying ST to the deep drilling data files from the available 40 of 64 CCD-quadrant sub-fields observing the galactic plane at [281, −24] deg in R.A. and decl. on 2018 July 13 as shown in Figure 1, we found 1168 objects including main belt asteroids in the core region of the asteroid belt, Trojan, Hilda, Hungaria, Phocaea asteroids, and near-Earth-asteroids (NEAs). Stacking up 133 30 s images enables us to detect signals as deep as about 23rd magnitude, while the typical detection threshold for each ZTF field is at about 20.5 magnitude with the ZTF-r filter (Masci et al. 2018). For this particular data set, the limiting magnitude at signal-to-noise ratio (S/N)=5 was 20.23 according to the fits metadata. We estimate about two-third of the objects are not detectable using a single frame. In the following sections, we will describe the ZTF data and data processing, present results, and conclude with discussions on the future of this approach.

Data and Method Description
ZTF focal plane is covered by 16 CCDs. Each CCD has 4 independent amplifier channels for each quadrant to output digital signals of array size 3080×3072 pixels. The plate scale is about  1 pixel −1 , so each quadrant covers a field of size about 0.73 deg 2 . The nominal exposure time is set to 30 s. The readout time is 8 s and the overhead is 2 s giving approximately a rate of 40 s per frame. The read noise is 10 e − per read. The detection limit per 30 s frame is about 20.5 mag for the ZTF-r band filter, limited by the sky background noise (Masci et al. 2018). The deep drilling data set that we analyzed was taken with the ZTF-r filter over the same field containing 133 frames giving total time span of about 5320 s. 8 The data files belong to the ZTF data product called Epochal-difference image files (Masci et al. 2018). For each CCD quadrant, a reference frame for the field was subtracted from each of the 30 s science frame based on the algorithm described in Zackay et al. (2016). The science frame is obtained from the raw frame after applying bias, flat field, and nonlinearity corrections. The streaks from aircraft and satellites, bright source halos, and ghosts are masked out. The reference frame provides a static representation of the sky and is generated by co-adding from 15 to 40 30 s images that pass specific data quality criteria. The frames we used have all the static objects in the sky subtracted down to the limiting magnitudes of the reference frames, which goes about 1.5 mag deeper than individual 30 s images from averaging 15-40 images (Masci et al. 2018).
As summarized in Figure 2, we start with the ZTF Epochaldifference image files and first use Matlab to prepare the data in the format suitable for our main ST search engine to run. We re-register the frames by shifting an integer amount of pixels according to the sky coordinate of the center of each CCDquadrant as recorded in the metadata. Shifting each frame removes the telescope tracking error, shown in Figure 3, so that all the 133 frames are aligned with respect to the sidereal. We then convert the ZTF data values into 16 bit unsigned integers as expected by our main ST search engine. The ZTF preprocessing subtracts out a static image of the sky for the field, so residual signals from static objects due to photon shot noise and atmospheric effects are left in the frames. This is different from our standard ST preprocessing where we mask out star signals typically down to S/N of 7 before subtracting static objects from each frame which leaves much less static object residuals in the frames. Because of the higher level of residual noises, to avoid generating excessive false positives, we do some extra preprocessing to regulate extreme pixel values. We set the following pixels to the sky background value: (1) pixels with counts greater than 600 DN (a value chosen to limit positive stellar residuals but not the signals from (2) pixels with counts that are 10σ below background or lower; and (3) pixels with integrated signal above S/N=10 after stacking up all the frames tracking sidereal. There is quite some room left for improving these ad hoc preprocessing to increase the detection sensitivity as discussed in Section 4.
The main ST search engine employs GPU-aided computation to perform shift-and-add of the frames in parallel for all the trial velocities in a 100×100 grid with grid spacing of 0.025 pixel/frame (~ 2. 3 hr −1 ) covering a range of velocity of ±0.75 deg day −1 in both R.A. and decl. Because half of the velocity grid spacing is the upper bound of tracking error for shift-and-add along each direction, we typically choose the spacing to be the rate at which the motion over the total integration time is 1-2 PSF size so the trailing loss is negligible. Our choice of 0.025 pixel/frame corresponds to a motion of about 1.5 PSF size per data cube. The PSF size is about  2. 3 limited by seeing. Our velocity range covers the speeds of the main belt asteroids at their opposition, which is between 0.13 and 0.33 deg day −1 . Trans-Neptune objects (TNOs) moves at rate in general less than  4 hr −1 or about  6 per data cube integration time. Our analysis of the deep drilling data set could have detected TNOs if they move more than 1 PSF size over the time of the data cube and are brighter than 23rd mag, but we did not find any in the field.
Using a single Nvidia P100 GPU, it takes about 220 s to process a 3080×3072 pixel CCD quadrant. For 64 quadrants, using 3 Nvidia P100 GPUs, we can process all the data within 4700 s, which is less than the total operational time of about 5320 s required for taking the deep drilling data, thus our shiftan-add implementation is viable for a quasi-real-time data processing. A Gaussian kernel with FWHM of 3 pixels is applied to improve the S/N. 9 For most of the CCD-quadrants, the detection threshold is S/N=10, where the noise level is determined by the spatial noise derived from spatial variations of pixel intensities in the frames, which varies between CCDquadrants. If the noises were uniform Gaussian, the expected false positive rate is close to zero. 10 However, a few CCDquadrants yield too many false positives using detection threshold of S/N=10 because noises are not uniform over the field and the software that we used to process the data only estimates an average spatial noise for each tracking velocity. For those data sets, higher S/N thresholds like 15 or 20 are used instead by manual adjustments. Signals above the detection threshold are found by the GPU-aided search and then clustered for post-processing.
In post-processing, each 133 frame data cube is broken into four sub-data cubes with equal number of frames. Therefore, each sub-data cube would have nominally 33 frames. For a few CCD-quadrants that do not have all the 133 frames (all have  9 Using a matched filter can improve this by about 5% but would require the detailed knowledge of the PSF shape. 10 For Gaussian noise, the probability to fluctuate above S/N is ( ) erfc S N 2 2, where erfc is the complementary error function. Therefore the expected false positive rate per data cube isŃ N N N x y Vx Vy x y is the array size of one CCD-quadrant and (N Vx , N Vy )=(100, 100) is the dimension of the velocity grid. For S/N=10, we get´´´3080 3072 100 100 ( )~´-erfc 10 2 2 7 10 13 . more than 100 frames), each sub-data cube would have slightly less frames. We require that the signal appears in at least 3 subdata cubes above S/N=5. When 4 sub-data cubes are combined, this gives S/N10). The PSF FWHM of each detected signal determined by fitting a Gaussian profile must also be within the range of [ ]   1. 2, 3. 5 or about 0.5-1.5 of average seeing limited PSF size (  2. 3). The field is quite crowded near the galactic plane with Galactic latitude −9°, so an asteroid has a good chance to get close to a bright star's halo, which is masked out as shown in Figure 5. Therefore, we only require the signal to appear in 3 sub-data cubes to include signals that are partially masked.
For ST search, we typically set a detection threshold much higher than the threshold used to detect signals in a single frame (e.g., S/N=5) to avoid false positives due to the lookelsewhere effect from integrating the frames for 100×100 different tracking velocities, i.e., there are 10,000 times more chance to get false positives ). For Gaussian noise, a detection threshold of S/N=10 would yield practically zero false positive rate; however, the dominant sources of noise are the residual stellar signals after subtracting the reference image. The noise is non-uniformly distributed, and the S/N is only calculated based on an average spatial noise level. The S/N noise quantity is calculated as the standard deviation of pixel-topixel variations over a CCD quadrant. If the track of the signal goes through a region that has excessive stellar residuals, the integrated signal could be biased high from including positive stellar residuals. We thus break the data cube into four sub-data cubes and require the signals to be in at least 3 sub-data cubes. Since the object moves, it would be unlikely that all the sub-data cubes are biased high by excessive stellar noises because the slowest object detected moves more than  6 per sub-data cube. Assuming two sub-data cubes are not significantly biased by excessive stellar noise fluctuation, the false positive rate is about 3080×3072×100×100×(erfc(5 2 )/2) 2 ∼0.0078 false detections per data cube. With the additional requirement of the appearance of the signal in the 3rd sub-data cube, we believe our false positive rate is negligible.

Results
In this section, we present examples of detected asteroid signals and statistics of the detections. The left plot in Figure 4 displays the signals in each sub-data cube with a color scale in units of temporal noise quantified as the median value of frame-to-frame variation of pixel intensities in four sub-data cubes when tracking the object at rate of [0.153, 0.373] pixel/ frame, or about [−35, −14]  hr −1 . This rate is typical for a main belt asteroid. This signal's magnitude is estimated to be 19.1. It is bright enough to be detected as a streak as shown in the right plot in Figure 4 even with the trailing loss when coadding all the 133 frames to simulate an integration of about 1.5 hr for ZTF tracking the sidereal. Although this can be detected as a streak (blue line) in the co-added images, for accurate astrometry, the whole data cube is needed to take the advantage of ST (Zhai et al. 2018). Similarly, Figure 5 displays another main belt asteroid of magnitude of 20.8, whose signals only show in three sub-data cubes because it ran into a masked region due to a bright star. Note that the signal in the second sub-data cube is also slightly lower because the signals of the asteroid in a few frames of the second sub-data cube are masked out. Even though the S/N of this object is about 50 in a sub-data cube, the streak of this object is hard to identify and is no longer significant due to the trailing loss in the right plot of Figure 5 when tracking the sidereal. Figure 6 displays in the same fashion an asteroid of magnitude of 22.4 whose signals consistently show in all four sub-data cubes. Note that the color scale is measured in units of temporal noise level, the actual S/N with respect to spatial noise level is less than 10. This discrepancy is because the signals are surrounded by relatively high residual stellar noises which fluctuate at level higher than the temporal noise. The signals along the track of the object (the blue line in the right image for tracking the sidereal) are completely buried in the noises. This shows the efficacy of ST. Figure 7 displays the histogram of the estimated magnitudes of the detected objects with a bin size of 0.5 stellar magnitude. We can see an increase of the population until 21st magnitude, showing the increase of the population of the asteroids with the decrease of sizes of asteroids. The drop of population in bins higher than 21st magnitude indicates an incompleteness of detection beyond 21st magnitude even though we can detect objects close to 23rd magnitude. The reason is because the limiting magnitudes varies a lot between CCD quadrants due to different spatial noise levels, PSF sizes, and data processing settings to avoid excessive false positives.   about 1.8 mag with the lowest quadrant limiting magnitude being 21.2, which is consistent with the fact that our detection is complete at magnitude of 21. The incompleteness comes also from the criteria of requiring signals to be in at least three subdata cubes because fainter signals tend to be more likely to become not significant if the signals are partially masked out or weakened by the fluctuations of spatial noises. Figure 9 displays histogram of the S/N of the detected objects where the noise level is measured by the standard deviation of the temporal noise. In the absence of stellar residuals and considering look-elsewhere effect for all the 64 CCD quadrants, using a detection threshold of S/N at 8 would give only a false positive rate of 3080×3072×100×100×erfc (8 2)×64≈0.0075. The decrease of population for S/N lower than about 30 suggests that we missed quite some objects with S/N between 8 and 30. Therefore, we expect lots of room to improve detection sensitivity to go beyond the current performance.
The sky rate of the detected objects are displayed in Figure 10. Different colors and dot sizes label different ranges of apparent magnitudes. The field center was at [281, −24] deg, observed on 2018 July 13, UTC 6-7, so observations were approximately along the anti-Sun direction with solar elongation about 170°. Assuming an object has a circular orbit with radius R o around the Sun in the ecliptic plane, viewed from the Earth at anti-Sun direction, the angular rate with respect to sidereal is ( ) deg/365.25 day, for R o /R e = 5.2 (Jupiter), 4 (Hilda), 2.5 (main belt core), we have sky rates of −0.13, −0.16, −0.24 deg day −1 respectively. Because we did not observe exactly at the anti-Sun direction and considering the small eccentricity of the Earth orbit, the actual rates are expected to be slightly slower than these values. The observation date is only a few weeks after the summer solstice, so the motion along ecliptic plane is very close to R.A. direction. Looking at the clusters in motion rate shown in Figure 10, we can easily identify asteroids of different families. The group near [−0.12, 0] deg day −1 are Trojan asteroids with with large decl. rates. They belong to the Hungaria family with semimajor axis between 1.78 and 2 au and with large inclination. At least one asteroid ((282212)=2001 XG 26 ) belongs to Phocaea family, which is between main belt asteroids and the Hungaria family with semimajor axis between 2.25 and 2.5 au and eccentricity greater than 0.1. In general, the sky rate of NEAs are not clustered and can have a wide range. There are at least two known NEAs (2018NE 4 , 2009 ) "re-discovered" with rates scattered outside the well-clustered groups. There should be many new asteroids detected here especially the fainter ones because optical survey of main belt asteroids is not complete beyond magnitude 21.5 (Ryan et al. 2015). Using the MPCChecker (https://minorplanetcenter.net/ cgi-bin/checkmp.cgi) tool provided by the Minor Planet Center (MPC), we found five asteroids with large negative decl. rates, marked by the two blue dots with decl. rate <−0.5 deg day −1 and three black dots with decl. rate <−0.2 deg day −1 , do not correspond to known asteroids. Some of them could be new NEAs.

Conclusions and Discussions
We have presented the data processing and results from applying the ST technique to a ZTF deep drilling data set to integrate over ∼1.5 hr. We found 1168 asteroids including main belt asteroids in the core region, Trojan, Hilda, Hungaria, and Phocaea family asteroids, as well as NEAs. These results demonstrate that this approach is powerful and productive for surveying all types of asteroids. The ST technique can significantly improve the sensitivity of a survey facility for detecting moving objects by combining multiple observations of the same field.
When a moving object is discovered, it is important to have its orbit determined so as to catalog it. This is usually done by the MPC, to which we have reported all the objects found here. MPC checks whether the observations are consistent with the orbits of known objects. For new objects (observations not attributed to any known object), MPC will make use of the observations to perform short term predictions of their ephemeris for facilities over the world to perform follow-up observations. However, because the deep drilling data set was not analyzed in time for follow-up observations, so our new objects only have observations over a time window of about 1.5 hr, thus are essentially lost, waiting for future recovery observations. Therefore, it is very important to process data sets and submit the new observations to MPC timely. Modern GPUs enable a quasi-real-time data processing using the ST, so this opens a new way for survey facilities to find asteroids, i.e., a deep drilling mode, by taking deep drilling data sets and making use of the ST, as has been advocated by Heinze et al. (2015). With many survey facilities are scanning the sky, same object tends to be found by multiple facilities. Going to deep drilling mode would be more efficient for finding new objects. Some of the objects discovered from deep drilling mode may to be too faint for most of the facilities to do follow-up, it may be useful to design the observation cadence so that the same field would be observed multiple times within a few weeks for a preliminary orbit determination. Note that ST enables to perform follow-up observations of multiple objects in a field conveniently because there is no need to track individual objects during observation. Another advantage of this approach is to produce more accurate astrometry (Zhai et al. 2018).
The work presented here is a preliminary study. More than one thousand objects per 1.5 hr is a very good rate of yield and worth pursuing further. In the future, we can search deeper by improving preprocessing to mask out more residual noises from the static objects. This can be done by using the raw ZTF data product and our standard data processing, which is fully automated. This would also make the manual adjustments of detection S/N thresholds for avoiding excessive false positives unnecessary. Note that majority of the data processing computation load is for performing the shift-and-add, so adding extra preprocessing would not increase processing time by more than 5%, thus not affect the capability of doing a quasi-realtime ST processing. Besides, our computation architecture allows adding more GPUs without requiring software change; so we are not limited by computation power yet.
There are tens of ZTF deep drilling fields with varying number of frames and field locations. We will work on one field that is less crowded with similar number of frames to allow us to study detection sensitivity when the signals of moving objects are less likely to be masked out due to star confusion. With improved preprocessing, we expect to see deeper beyond 23rd magnitude and significantly increase completeness in magnitude range 21-23 for surveying slowly moving objects. The majority of these fainter objects to be detected will be very likely new detections.