Tails: Chasing Comets with the Zwicky Transient Facility and Deep Learning

The ZTF ¹² is a state-of-the-art robotic time-domain sky survey capable of visiting the entire visible sky north of −30° decl. every night. ZTF observes the sky in the g, r, and i bands at different cadences depending on the scientific program and sky region (Bellm et al. 2019a; Graham et al. 2019). The 576 megapixel camera with a 47 deg² field of view, installed on the Samuel Oschin 48 inch (1.2 m) Schmidt Telescope, can scan more than 3750 deg² per hour, to a 5σ detection limit of 20.7 mag in the r band with a 30 s exposure during a new moon (Masci et al. 2019; Dekany et al. 2020).

The ZTF partnership has been running a specialized survey, the Twilight Survey (ZTF-TS) that operates at solar elongations down to 35 degrees with an r-band limiting magnitude of 19.5 (Ye et al. 2020; Bellm et al. 2019b). ZTF-TS has so far resulted in the discovery of a number of Atira asteroids (orbits interior to Earth's) as well as the first inner-Venus object, 2020 AV2 (Ip et al. 2020). Motivated by the success, ZTF-TS will be expanded in Phase II of the project, which commenced in 2020 December.

Comets become more easily detectable when close to the Sun as they become brighter and start exhibiting more pronounced coma and tails. Furthermore, it has been shown that the most detectable direction of approach of an interstellar object is from directly behind the Sun because of observational selection effects (Jedicke et al. 2016) and the fact that this direction has a greater cross section for asteroids to bend around and pass into the visibility volume (Engelhardt et al. 2017; Do et al. 2018).

Tails automates the search for comets with detectable morphology. While trained and evaluated on a large corpus of ZTF data, in this work we focus on Tails' performance when applied to the ZTF-TS data.

To build a seed sample for labeling, we first identified all potential observations of known comets conducted with ZTF from 2018 March 5 to 2020 March 4, based on their predicted position and brightness. The code for accomplishing that is based on the Python libraries pypride (Duev et al. 2016) and solarsyslib (Jensen-Clem et al. 2018) and uses the comet ephemerides obtained from the Minor Planet Center (MPC) ¹³ for a coarse search, followed by a JPL Horizons ¹⁴ (Giorgini et al. 1996) query for precision.

To provide more contextual information, epochal image data are supplemented by properly aligned reference images of the corresponding patches of sky and difference (epochal minus reference) images generated with the ZOGY algorithm (Zackay et al. 2016), all produced by the ZTF Science Data System at Caltech's IPAC (Masci et al. 2019). Finally, we generate image triplet cutouts 256 by 256 pixel in size, which in angular measure translates into 43 by 43 at ZTF's pixel scale of 101/pixel.

We selected over 60,000 individual observations with the total comet magnitude ranging from 10 to 23 (as reported by JPL Horizons; see Figure 1), out of which about 20,000 were sourced for manual annotation. This resulted in an initial sample of 3000 examples with identifiable morphology.

Zoom In Zoom Out Reset image size

We also compiled a set of approximately 20,000 negative examples consisting of point-like cometary detections, patches of sky with no identified transient or variable sources, CCD-edge cases, and a wide range of real (point-source) transient and bogus (e.g., artifacts due to bright stars, optical ghosts, and "dementors") samples from the Braai data set (Duev et al. 2019).

To expand the data set, we then assembled a standard ResNet-based (He et al. 2015) classifier for comet identification. With this basic classification model, we ran several rounds of an active-learning-like procedure, where we would first train the classifier, evaluate it on the whole data set, sample both confident predictions and the cases close to the classifier's decision boundary, manually inspect and label those examples, and add them to the training set. Roughly 2000 positive and 2000 negative examples were added to the training set via this method.

The resulting training data set contains about 5000 positive and 22,000 negative examples (see Figure 2). Each triplet in the set has been assigned a label [p_c , x, y], where p_c marks the presence of a comet in the image and x, y ∈ [0, 1] is the relative positions of the comet's center of mass, as reported by JPL Horizons. For positive examples this translates into [1, x_JPL, y_JPL], and for negative ones [0, ?, ?], where the question marks denote that these do not affect the loss in this case.

Zoom In Zoom Out Reset image size

Tails adopts a custom architecture (see Figure 3) based on EfficientDet D0 (Tan et al. 2019), a variant of a state-of-the-art architecture designed for object detection—a computer-vision technique for the identification and location of objects in image data.

Zoom In Zoom Out Reset image size

This architecture delivers best-in-class object detection efficiency and performance across a wide range of resource constraints. This is achieved by using EfficientNet—state-of-the art backbone networks for feature extraction, a weighted bidirectional feature pyramid network (BiFPN), which allows easy and fast multiscale feature fusion, and a compound scaling method that simultaneously and uniformly scales the resolution, depth, and width for all backbone, feature, and location/class prediction networks (Tan et al. 2019).

The use of a BiFPN, which effectively represents and processes multiscale features, makes this architecture particularly well suited for the problem of morphology-based comet identification and localization.

A batch of triplet image stacks (n_b , 256, 256, 3) in size, where n_b is the number of stacks in the batch, is passed through an EfficientNet B0 backbone (Tan & Le 2019). The extracted features from the last five blocks/levels of the network are passed through the BiFPN. The resulting five output tensors denoted by the colored circles in Figure 3 are fed into the head network, which outputs the probability of the image containing a comet p_c and its centroid's predicted relative (x, y) position. ¹⁵

We defined the loss function as

$\begin{eqnarray}&&L={w}_{c}\cdot {L}_{c}+{w}_{p}\cdot {L}_{p},\end{eqnarray} \tag{ 1 }$

where L_c denotes the binary cross-entropy function for the label c (1—there is a comet in the image, 0—there is no comet) and the predicted probability p_c . If ⌊p_c ⌉ = 1, L_p is computed as an L₁ loss for the relative position (x, y) and its prediction (x_p , y_p ) with a small L₂ regularizing term (with = 10⁻³), and w_c and w_p denote the weights of the two terms, respectively,

$\begin{eqnarray}\begin{array}{rcl}{L}_{c} & = & \displaystyle \sum c\cdot \mathrm{log}({p}_{c})+(1-c)\cdot \mathrm{log}(1-{p}_{c})\\ {L}_{p} & = & \displaystyle \sum \lfloor {p}_{c}\rceil \cdot \left(| x-{x}_{p}| +| y-{y}_{p}| \right.\\ & & \left.+\,\epsilon \cdot \sqrt{| x-{x}_{p}{| }^{2}+| y-{y}_{p}{| }^{2}}\right).\end{array}\end{eqnarray} \tag{ 2 }$

We employed the Adam optimizer (Kingma & Ba 2014), a batch size of 32, and a 81%/9%/10% training/validation/test data split. For data augmentation, we applied random horizontal and vertical flips of the input data; no random rotations and translations were added. We note that the test/validation sets did not contain augmented data from the training set. We used standard techniques to maximize training performance: if no improvement in validation loss was observed for 10 epochs, the learning rate was reduced by a factor of 2, and training was stopped early if no improvement was observed for 30 epochs.

The EfficientNet's weights were randomly initialized. ¹⁶ We first set w_c = 10, w_p = 1 to allow for a fast convergence of the feature-extracting part of the network. To fine-tune the performance, we trained Tails on a balanced data set setting w_c = 1.1, w_p = 1 and monitored the validation loss for early stopping, then bumped w_p = 2 and monitored the validation positional rms error, and finally, we added the omitted negative examples and again monitored the validation loss for early stopping.

The resulting classifiers were put through the same active-learning-like procedure as was employed in the initial data set assembly, using several months of ZTF-TS data.

On 2020 October 7, Tails discovered a candidate that was posted to MPC's Possible Comet Confirmation Page (PCCP) ²⁰ as ZTFDD01 (see Figure 7). It was later confirmed to be a long-period comet and designated C/2020 T2 (Palomar), marking the first DL-assisted comet discovery (Duev et al. 2020). The candidate was found in the Twilight survey data; it was at 19.3 mag in the ZTF r band. The FWHM of the object was approximately 25–3'', compared to nearby background stars that have an FWHM of ∼2''. The object showed a tail extending up to 5'' in the westward direction. Table 1 summarizes the orbital elements of C/2020 T2 provided by the MPC and Figure 8 shows its orbit as of the discovery date.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine if Tails could have discovered C/2020 T2 before 2020 October 7, we searched the ZTF archive for all Twilight Survey data covering the ephemeris position of the comet with the ZChecker software (Kelley et al. 2019). Eleven nights of data were found between 2020 June 11–20 (evening twilight) and October 7–21 (morning twilight). The comet was in conjunction with the Sun between the two sets, and not observable by ZTF. We measured the brightness of the coma in 4 pixel radius apertures, and aperture corrected the photometry according to the ZTF pipeline documentation. The data are shown in Figure 9. Typical seeing was 2'' in June, and the comet was very faint (r = 20.2 mag), near the single-image detection limit (r = 20.4–20.9 mag, 5σ point source), and had no morphological features for Tails to pick up. Thus 2020 October 7 was really the first opportunity for Tails to discover the comet.

Zoom In Zoom Out Reset image size

A comet candidate was identified by a combination of Tails and the ZTF Moving Object Detection Engine (Masci et al. 2019) on 2020 January 4 UTC and was submitted to the PCCP as ZTF0Ion (see Figure 10). It was later identified as a recovery of comet P/2016 J3 (STEREO) and given the designation P/2016 J3 = P/2021 A3 (STEREO; Bolin 2021). P/2021 A3 was identified in the evening Twilight survey data at r 19.3 mag with a clearly extended appearance scoring 0.9 with a coma ∼10'' wide and a tail extending past 20'' in the northeast direction.

Zoom In Zoom Out Reset image size