2.1. SkyMapper telescope

SkyMapper is a 1.35 m modified Cassegrain telescope located at Siding Spring Observatory (SSO) in New South Wales, Australia (Wolf et al. Reference Wolf2018). The telescope has a wide field of view of $2.34 \times 2.37\, \mathrm{deg}^2$ , a uvgriz filter set (Bessell et al. Reference Bessell, Bloxham, Schmidt, Keller, Tisserand and Francis2011) and a mosaic of 32 $2\mathrm{k}\times4\mathrm{k}$ CCD detectors with a pixel scale of $0.5 \mathrm{arcsec pixel}^{-1}$ . It is owned and operated by the Australian National University (ANU) and, most importantly for EM follow-up, it is a robotic facility.

Figure 1. Top left: LIGO/Virgo probability sky map for S190814bv produced with the ligo.skymap module (Singer et al. Reference Singer2016). The inset shows the most probable area for the optical counterpart. Darker colours correspond to higher probability sky regions, and contours enclosing regions with 50%, 90%, and 99% of the probability are indicated. Bottom left: probability map convolved with the coverage of reference images in SkyMapper DR2. Bottom right: zoomed-in map of the 20 highest-probability fields selected for the search; here, one field alone has $\sim 39\%$ probability of containing the GW source. Symbols are bright galaxies from the 2MASS redshift survey. Top right: Observability plot for the top 20 fields with telescope altitude, night-time range and Moon separation.

The main purpose of the telescope is the SkyMapper Southern Survey (SMSS), a hemispheric sky atlas which has been underway since 2014 (DR1: Wolf et al. Reference Wolf2018; DR2: Onken et al. Reference Onken2019). For point sources with $\mathrm{SNR}>5$ , the expected photometric depth of single-epoch 100-s images is u = 19.5, v = 19.5, g = 21, r = 20.5, i = 20, and z = 19 mag. Alongside, SkyMapper has been used to search for extragalactic transients, including low-redshift Type Ia supernovae (Scalzo et al. Reference Scalzo2017; Möller et al. Reference Griffin2019) as well as optical counterparts to GW events (GW170817: Abbott et al. Reference Abbott2017b; Andreoni et al. Reference Andreoni2017) and fast radio bursts (e.g. Petroff et al. Reference Petroff2015; Farah et al. Reference Farah2018; Reference Farah2019; Price et al. Reference Price2019; Chang et al. Reference Chang, Wolf, Onken, Luvaul, Gupta and Flynn2019a). Different needs made the SMSS and the transient searches develop separate data reduction software: the Science Data Pipeline (SDP: Luvaul et al. Reference Luvaul, Onken, Wolf, Smillie, Sebo, Lorente, Shortridge and Wayth2017; Wolf et al. Reference Wolf2018) for the SMSS and the Subtraction Pipeline (SUBPIPE; Scalzo et al. Reference Scalzo2017) for the Transient Survey. SUBPIPE was used for EM data analysis before before LIGO/Virgo ER14 (Engineering Run 14), but for the O3 run we updated it to an automated real-time pipeline (AlertSDP: Alert Science Data Pipeline) that borrows many features from the SDP (see Sections 2.2 and 2.3 for details).

2.2. Alert response

We continuously listen to the live stream of LIGO/Virgo public alerts for compact binary merger candidates.Footnote a The stream is distributed through the Gamma-ray Coordinates Network (GCN) using the pygcn Python module.Footnote b We developed a robotic alert handler that extracts relevant information, ingests the GW event into our database, downloads the HEALPix 3D localisation map (skymap), prioritises areas for follow-up, and generates a list of observations for SkyMapper.

The first preliminary GCN notice is automatically issued for a superevent within 1–10 min after the GW trigger. From this notice, our robotic handler initiates rapid-response search observations by convolving the probability in the BAYESTAR skymap (Singer & Price Reference Singer and Price2016) with the tile pattern of the SkyMapper Southern Survey (Onken et al. Reference Onken2019), ranking the fields by probability integrated per field and then selecting the top 20 fields with existing i-band reference images ( $\sim100\, \mathrm{deg}^2$ total area). Figure 1 shows the sky maps for a real alert and an observability map of the selected target fields.

Figure 2. The i-band light curve for the GW170817 kilonova, at different distances: at the true distance of 40 Mpc (top), shifted to 100 Mpc (middle) and 200 Mpc (bottom); solid lines are power law decay fits. The dashed line at $i_{\textrm{AB}}=20$ marks our typical 5- $\sigma$ magnitude limit in 100 s exposures (data were taken from various literature sources; AST3-2: Hu et al. Reference Hu2017 ; B&C: Utsumi et al. Reference Utsumi2017; DECam: Cowperthwaite et al. Reference Cowperthwaite2017; Gemini: Kasliwal et al. Reference Kasliwal2017; LaSilla: Smartt et al. Reference Smartt2017; LCO: Arcavi et al. Reference Arcavi2017; Magellan: Shappee et al. Reference Shappee2017; Pan-STARRS: Smartt et al. Reference Smartt2017; REM: Pian et al. Reference Pian2017; SkyMapper: Andreoni et al. Reference Andreoni2017; Swope: Drout et al. Reference Drout2017;T80S: D az et al. 2017; VLT: Tanvir et al. Reference Tanvir2017; VST: Pian et al. Reference Pian2017).

In this search stage, we obtain two $100\,\mathrm{s}$ exposures in i-band for each field, separated by $\sim 8\,\mathrm{min}$ . Requiring a twin detection in the two images eliminates moving objects from the candidate list. The typical rate of motion for main-belt asteroids, 0.625 arcsec $\mathrm{min}^{-1}$ near opposition (Jedicke Reference Jedicke1996), provides a sufficient $\sim 5$ arcsec spacing between observations. With a depth of $i\approx 20\,\mathrm{mag}$ , we can identify kilonovae that are $\sim 3\,\mathrm{mag} (16\times$ ) fainter than the kilonova of GW170817 was 10 h after the GW trigger. Hence, we could detect kilonovae for GW170817-like events out to $4\times$ the distance of GW170817, or 160 Mpc. Figure 2 shows the lightcurve of the GW170817 kilonova shifted to different distances up to 200 Mpc. The declining nature of the i-band light curve suggests that kilonovae may be more luminous during the first 10 h. We thus assume that we might be able to detect early kilonova emission, just hours after a BNS merger, out to distances of 200 Mpc or beyond.

We search in i-band for three reasons: (i) kilonovae are expected to be red due to the large opacity of neutron-rich ejecta, although at early times (within an hour of the merger) blue emission from hotter, polar ejecta may be expected (e.g. Troja et al. Reference Troja2018), (ii) the i-band covers the largest sky area with deep reference images in the SMSS that can be used for real-time image subtraction, and therefore provides the greatest search volume for kilonovae, and (iii) the dominant source of transient events, flares on M stars too faint to be seen in quiescence, is nearly irrelevant in i-band (Chang, Wolf, & Onken Reference Chang, Wolf and Onken2020). The aim of the search phase is to find a counterpart as fast as possible, report it to facilities around the world and get SkyMapper itself into monitoring mode. We expect that the search phase provides a full set of possible counterparts within 4 h (limited by image processing) after the GW trigger, if the trigger occurs during darkness and the sky is clear. Later, possibly within 4 h for BNS or NSBH events, an updated LALInference skymap (Veitch et al. Reference Veitch2015) will be distributed, including an updated sky localisation and source classification. This sometimes leads to a significant change in the sky localisation, and the SkyMapper observing plan may be modified and re-executed as a result.

At any time, once a position of a likely kilonova transient is identified, we can manually switch from the search phase to a continuous high-cadence monitoring of the new source. Observations of an early kilonova with a cadence of 2 min would reveal structure in the lightcurve arising from shocks induced by the kilonova ejecta (e.g. Metzger Reference Metzger2019). We will also get high-cadence, multi-band light curves for several consecutive nights after any kilonova discovery. The monitoring strategy will change depending on the colour, luminosity, and fading timescale of the kilonova candidates. If the optical counterpart gets identified by other groups before SkyMapper can observe, we will do only light-curve monitoring, as was the case for GW170817 (Andreoni et al. Reference Andreoni2017).

2.3. AlertSDP: Alert science data pipeline

The AlertSDP is a new pipeline created from two existing software packages, which are the SkyMapper Science Data Pipeline, SDP, and SUBPIPE, which was used for discovering low-redshift supernovae (see Section 2.1). The SDP has provided the general process control framework as well as improved calibration and masking treatment, while the subtraction and transient classification have evolved from SUBPIPE components. First, each raw image is preprocessed and calibrated. After that, i-band reference images for image subtraction are selected, and transient candidates detected on difference images are uploaded into a database. The database also includes external catalogues to assist transient classification.

2.3.2. Selection of reference images

Reference images are required for image subtraction and real-time transient searching. We take advantage of the nearly all-Southern-sky coverage ( $>\!98 \%$ ) of the deeper Main Survey exposures from the SMSS in i-band ( $\lambda_\mathrm{centre}/\Delta \lambda = 779/140 \mathrm{nm}$ ). Figure 3 shows the area of sky covered by SkyMapper DR2 (Onken et al. Reference Onken2019) and by DR3, which was released to the Australian community in February 2020 in the late stages of the O3 run. DR3 expands the deep coverage by including survey observations from 2014 March to 2019 October.Footnote c

Figure 3. SkyMapper i band coverage in DR2 (top) and DR3 (bottom). Darker colours resemble deeper reference images used for subtraction.

From any new image, we subtract each available reference image. The best possible reference image combines a small PSF and a large overlap with the new image. The reference images come with a range of point-spread functions (PSF), and because the SMSS obtains repeat images with dithering to cover the sky homogeneously, the best-available reference PSF changes discontinuously with sky location. Small overlaps lead to badly determined convolution kernels and thus bad subtractions. Hence, we require at least 15% overlap area for mosaic frames and at least 5% overlap for each individual CCD. Following a common strategy of other transient surveys, we require that the reference image should be taken at least two weeks prior to a given new epoch.

2.3.3. Subtraction of images

We perform image subtraction on each field of overlapping new/reference image pairs. First, we resample each of the reference mosaic frames onto the new image with SWarp (Bertin et al. Reference Bertin, Mellier, Radovich, Missonnier, Didelon, Morin, Bohlender, Durand and Handley2002). Next, we convolve those pairs of images to a common PSF using HOTPANTS (Becker Reference Becker2015), which implemented the popular image subtraction algorithm by Alard & Lupton (Reference Alard and Lupton1998). Solving for the convolution kernel is a crucial step to equalising the PSFs of the reference and new images. The position-dependent PSF variation in both images is modelled as a linear combination of basis functions, and by default we choose a 2D polynomial of order two. Next, the flux level of the subtracted image is normalised to that of the reference image. Since the photometric calibration is directly tied to the SMSS photometry, this approach has the advantage of allowing explicit calculation of zero-point corrected magnitudes. The preliminary calibrated magnitude and error used for lightcurve construction is

\begin{align} m_{\rm sub} &= -2.5 \log f_{\rm sub} + {\rm ZP_{\rm ref}}, \notag \\[4pt] \delta m_{\rm sub} &= \sqrt{(1.0857 \times \delta f_{\rm sub}/f_{\rm sub} )^{2} + {\rm \delta ZP_{\rm ref}}^{2}}, \notag\end{align}

where ${f_{\rm sub}}$ is the difference flux on the subtracted image and ${\rm ZP_{\rm ref}}$ is the zeropoint of the reference image. Here, we use the 15 arcsec (30-pixel-diameter) aperture as a total-magnitude reference. The zeropoint error $\delta {\rm ZP_{\rm ref}}$ contributes little, since only images with robust zeropoints are included in SMSS data releases.

Finally, we run SExtractor (Bertin & Arnouts Reference Bertin and Arnouts1996) on the subtracted images to produce a list of transient candidates and their associated metadata, using a low detection threshold of $1.5\sigma$ above the background. Gaussian filters in SExtractor are applied to an image prior to the detection of sources. Only sources with a signal-to-noise ratio below 2 are referred to as sub-threshold detections.

2.4. Transient classification

Until the end of O3, we classified transient candidates with a random forest (RF) classifier trained on earlier supernova survey data (Scalzo et al. Reference Scalzo2017). We computed a set of features derived from difference images of individual candidates, similar to those proposed in Bloom et al. (Reference Bloom2012). The main difference was that three additional checks were made on each candidate: (i) we removed any detections where the measurement of shape parameters (e.g. FWHM, elongation) was significantly larger than the median quantities for sources found in the new science images, and (ii) we matched each of them to the nearest object ( $<30 \mathrm{arcsec}$ ) seen on the reference frame, which will be associated with host galaxies or bright stars that might be poorly subtracted. Then, the RF model assigned a real-bogus score (RBscore) from 0 (artefacts = bogus) to 100 (transients = real) to all the detected candidates.

In this work, we added a new classifier (XGBoost) and new training sets to exploit our improved image processing. A combination of classifiers is often more accurate than a single classifier, and thus we introduce a new metric (Tscore) that combines RBscore and XGBoost score to compensate for the weakness of each classifier (see Section 3 for the full description and performance).

For each transient candidate, we identify whether it may be a known variable or Milky Way source by cross-matching it against catalogues of stars, variable objects and quasars, using a 3 arcsec search radius. We also reject known solar system objects with the SkyBoT cone-search service (Berthier et al. Reference Berthier, Vachier, Thuillot, Fernique, Ochsenbein, Genova, Lainey, Arlot, Gabriel, Arviset, Ponz and Enrique2006). Our preselected categories are as follows:

1. • We define a ‘Star’ sample using parallaxes and proper motions (PPM) from Gaia eDR3 (Gaia Collaboration et al. Reference Gaia2016, Reference Gaia, Brown, Vallenari, Prusti, de Bruijne, Babusiaux and Biermann2020), requiring that the significance of a PPM signalFootnote d is $S_{\rm PPM}>3$ .

2. • For Gaia eDR3 objects with lower PPM SN or no PPM information, we use a ‘GaiaSource’ label.

3. • We use the label ‘Var’ for all sources cross-matched to the AAVSO International Variable Star Index (Watson, Henden, & Price Reference Watson, Henden and Price2006; version 2020 October 19). It contains mainly known variable stars and a few QSOs.

4. • We assign a label of ‘Quasar’ to sources matched to the Milliquas catalogue (Flesch Reference Flesch2015, version 7.0a).

5. • For detections that are too close to very bright stars and carry a risk of being an optical reflection (defined by the algorithm in item 3 of Section 2.6.1 from Onken et al. Reference Onken2019), we use a ‘BrightStar’ label.

In addition, image mask data allows us to reject sources as being likely spurious, when they could be affected by bad pixels, neighbouring saturated pixels, or cosmic ray hits. Only candidates with Tscore greater than 30 (see Section 3) are considered plausible candidates and presented for visual inspection. If the number of plausible sources in a single CCD image exceeds 50, we reject all detections in that CCD because the quality of subtracted image is likely to be poor.

2.5. Further follow-up mechanisms

We manage data sets of transient candidates, building on a web user interface developed by Scalzo et al. (Reference Scalzo2017) based on the Django framework. Following the initial automated typing of the candidates, the next step requires a human to visually inspect and classify the remaining candidates in a web interface. All plausible candidates are assigned a unique name. Since we visit any field at least twice during the search, we prioritise candidates in the visual examination that are detected at least twice. We keep a list of ‘active’ candidates, which human vetting has identified as relevant for spectroscopic or photometric follow-up. If another facility reports a kilonova discovery, we can add such an object manually.

The web interface displays light curves, thumbnails of new, reference and subtracted images, and information from external services. Because of the rapid nature of the kilonova evolution and the robotic nature of SkyMapper operation, we added features to the web interface that trigger interaction with the telescope scheduler for follow-up modes that we anticipate to use during the next GW observing run (O4):

1. • Candidate Dweller: this feature can be used to monitor one or several objects in one of three possible modes: (1) The ‘one-off’ mode triggers just three exposures in the filter sequence i-u-i using an exposure sequence of 100-300-100 s (default but changeable). This mode is designed to probe the colour of a source quickly after an initial detection is made to help assess its likelihood of being a kilonova from a BNS merger. Two further modes with different sampling patterns can be used for continuous monitoring. (2) The ‘sampling’ mode probes the lightcurve and its colour evolution with alternating 100 s i band and 300 s u band exposures. And (3) the ‘intense sampling’ mode takes a series of ten consecutive 100 s i-band images followed by one 300 s u-band image, which provides the highest possible cadence, while recording colour evolution at lower cadence. Both sampling modes will continue observing the chosen target(s) until it sets or the dome closes.

2. • Alert Schedule Cleaner: this feature simply removes all observations in the queue, aborting any previously committed sequence and allowing a fresh start.

3. • AlertSDP Status: this web page reports the execution status of individual pipeline tasks (Section 2.3) and includes execution times for all jobs in the workflow.

We also trigger prompt Target-of-Opportunity spectroscopy at the ANU 2.3m Telescope to verify the physical nature of relevant candidates and inform further SkyMapper activities.

3.1. New XGBoost classifier

We adopt XGBoostFootnote e using an ensemble of decision trees. It uses a gradient boosting algorithm (Friedman Reference Friedman2000) to minimise a loss function when adding new decision trees. Unlike random forest methods, which train each tree independently, gradient boosted trees are built sequentially such that each subsequent tree aims to reduce the errors from its predecessors. The accuracy of classification is improved as more trees are added to the model, although a large complexity of the trees can lead to overfitting. However, XGBoost provides additional regularisation hyperparameters that can help reduce model complexity and guard against overfitting. Therefore, we use a column subsampling option to ensure that it uses a random subsample of the training data prior to growing trees.

The robustness of the classifier model is mainly limited by the available training data of candidates with known class label (e.g. Brink et al. Reference Brink, Richards, Poznanski, Bloom, Rice, Negahban and Wainwright2013; Scalzo et al. Reference Scalzo2017). The labelled data from the earlier supernova survey was unfortunately not useful in this regard, because it was predominantly obtained in bad seeing conditions. A new training set was generated by randomly selecting 1 000 SMSS DR3 images in i-band, excluding low galactic latitudes of $\mid b\mid <15^{\circ}$ . These data represent a large range of observational conditions and image quality from the survey. In this data set, we searched for transient candidates and used the original RF classifier to filter the list of millions of candidates. We then eyeballed candidates with RBscore > 40, and labelled them as real or bogus, where the latter category includes bad subtractions, cosmic ray hits and warm pixels. We also added to the real set known asteroids, variable stars, quasars and a small number of candidates projected onto galaxies from the 6dFGS and 2MASS XSC catalogues, as these may appear like typical host galaxies of kilonovae. Finally, we added to the bogus set a random subset of candidates with RBscore < 30, provided they were not associated with known objects of those types mentioned above that may be genuinely variable; this sample mostly includes subtraction artefacts such as residual features around bright stars, bad pixels and cosmic rays. The real-labelled class has much fewer instances than the bogus-labelled class, with a number ratio close to 1:10.

Table 1. Explored and chosen (bold) XGBoost Hyperparameters. Note that the results did not vary strongly with parameter changes

For training the XGBoost classifier with the similar input features used in Scalzo et al. (Reference Scalzo2017), we split the sample into training and testing sets of 43 593 and 71 307 candidates, respectively. One issue is that our binary classification does not have a balanced number of instances in the training set. This requires the use of a stratified sampling strategy to learn the features of each class equally. Next, the classifier has a list of hyperparameters that require fine-tuning in order to derive the best-possible model. We select test values of hyperparameters from a set of points in a coarse grid (see Table 1). We focused on hyperparameters that tend to have a high impact on the classification, such as control overfitting, learning rate, and complexity of the trees. Briefly the parameters are:

1. • learning_rate: step size shrinkage used in update to prevent over-fitting,

2. • max_depth: maximum depth of a tree,

3. • n_estimators: the number of trees in our ensemble,

4. • colsample_bytree: the subsample ratio of columns when constructing each tree,

5. • lambda: L2 regularisation term on weights,

6. • alpha: L1 regularisation term on weights.

The final parameters were chosen as the model with the lowest classification error, i.e. with both high completeness and high purity (see bold figures in Table 1).

3.2. Validation sample

To enable an accurate and unbiased assessment of classifier purity (Section 3.3), we need to define a ‘validation sample’—a random sample of objects with the weight factors based on the observed magnitude and RB scores. In the general case of an unweighted random sample, high RB scores are less common than low ones. Also, it is required to keep reasonable statistics for the rare, bright objects. We note that the distribution of classes in the validation set is unbalanced and neither reflects those in the training set perfectly nor those expected during the real transient search. The role of the validation set is to lead the parameters of the classifier towards best performance on real data. Hence, the distribution of classes in the validation set should ideally reflect that of the classes in the test set, so that the performance metrics will be similar on both sets. In other words, the validation set should reflect the expected data imbalance. The imbalance in our validation set thus leads to suboptimal performance. In Section 4 we test its performance in a real transient search during the O3 run.

From a set of 149 DR3 i-band images, we obtained 11 790 sources with RBscore larger than 30 without any further filtering, and we selected 4 500 candidates for further visual inspection in a manner that sampled the range of source parameters. We label the sources irrespective of whether they have an automatic label or not. Table 2 summarises the steps we used to select pure transient candidates by removing known image artefacts (cleaned sample 1) and preselected categories (cleaned sample 2). Sample 1 can evaluate basic system performance, while sample 2 represents transient candidates that remain unexplained after automatic association with a known variable source and need to be presented for human vetting when hunting for real kilonovae.

Table 2. Sample selection criteria for purity test

Our main interest is in classifying either isolated (with no apparent host galaxy) or supernova-like transients with high completeness. To measure the classifier completeness (see Section 3.4), we initially used a total of 5 194 asteroid detections that were identified by SkyBot and 443 detections for 26 supernovae (16 SN Ia, 7 SN II, 1 SN IIn, 1 SN Ibc, and 1 SN Ic) that had been followed-up or discovered by the SkyMapper Transient Survey (Scalzo et al. Reference Scalzo2017; Möller et al. 2019).

We additionally use the Open Supernova Catalogue (Guillochon et al. Reference Guillochon, Parrent, Kelley and Margutti2017) to collect spectroscopically confirmed SNe with distances less than 250 Mpc. In order to match known SNe against DR3, we take a broad date range of $\pm30 \mathrm{d}$ from the SN discovery date as our reference epoch. With this selection cut, there are 429 i-band images of 318 SNe, where the position and epoch were matched with sources in the SMSS DR3 catalogue within 5 arcsec. Some of those DR3 sources are the host galaxies of SNe, or otherwise chance superpositions of unrelated sources. After running the image set through the AlertSDP pipeline, we recover 153 detections for 117 SNe which lie in the magnitude range between 15.5 and 20.5 in i band. Since the SMSS is focused on covering the entirety of the Southern Sky rather than frequently repeating a given area, most of the SNe only appear in a single exposure. By type, this sample contains 69 Type I SN (64 Ia, 2 Ib, 2 Ic, 1 Ib/Ic), 45 Type II SN (36 II, 4 IIP, 2 IIb, 3 IIn), and 3 unclassified ones, in a variety of host galaxies. Thus, our completeness test with the Open Supernova Catalog sample will be less affected by host galaxy selection effects than that of the SMT SNe sample. Figure 4 shows thumbnail images of representative SN examples in different redshift ranges. We assume that this validation sample contains a more representative sample of contaminants, but also represents the outcome of the pipeline more accurately than the training data used in the model development stage.

Figure 4. Thumbnail images of two supernovae in the SMSS DR3 i-band data set, showing the new, reference, and subtraction image from left to right (size: $1\, \textrm{arcmin} \times 1\, \textrm{arcmin}$ ). Top: Type II SN 2019ejj at $d=13\, \textrm{Mpc}$ , i = 17.1. Bottom: Type Ia-91T SN 2019ur at $d=250\, \textrm{Mpc}$ , i = 18.8.

Figure 5. Purity vs completeness curves for the two cleaned samples described in Section 3.2. We compare the performance of three ensemble scores: RBscore (grey filled circle), XGBscore (grey open circle), and our new metric, Tscore (black filled triangle). The text in the panels refers to purity (P) and completeness (C) scores at a threshold of $t=30$ , which we adopt for our transient search.

Figure 6. Completeness test with asteroid (left), SMT SN (middle), and DR3 SN (right) samples as a function of magnitude. We compare the three different ensemble scores: RBscore (grey filled circle), XGBscore (grey open circle), and Tscore (black filled triangle). The same axis ranges are used in each panel.

3.3. Purity test

As our two classifiers have different strengths, a combination of them can lead to better performance. XGBoost is much better at recognising bad subtractions with negative pixel values, whereas RF works better at recognising that CR hits and warm pixels are not real objects. While this may be the result of different training, we define a new metric here, Tscore, that combines the two machine learning models (known as an ensemble of classifiers; see Chapter 34 in Dietterich Reference Dietterich2000 for instance), using the simple rule:

\[ {\rm{Tscore}} = \frac{\rm RBscore + \rm XGBscore}{2},\]

which gives both classifiers a similar weight.

For any threshold t, we define classifier purity:

\[ \mathrm{Purity}= \frac{\rm TP}{\rm TP + \rm FP},\]

where TP is the number of true transients with any scores greater than t, and FP is the number of false positives that were incorrectly considered as real. Varying the threshold of a binary classifier usually trades off better completeness against better purity. In Figure 5, we show curves of completeness vs purity in our cleaned samples as a function of threshold and compare our two classifiers as well as the combined Tscore. In cleaned sample 1, XGBscore is an improvement over RBscore at every completeness $>90\%$ ; for further applications, we choose a threshold of $t=30$ , which delivers 92% purity and 97% completeness. However, at higher purity the originally used RBscore appears more complete. Our new metric then combines the advantages of both and is in all parts of the curve at least as good as either the RB or XGB classifier. In cleaned sample 2 (right panel of Figure 5) all transient candidates explained by known objects have been removed and only those in need of human inspection are left over; here, a threshold of $t=30$ is still a good choice and provides a purity above 80% with 99% completeness. Based on this, we choose the ensemble-based Tscore classifier with a threshold of $t=30$ as our final classifier. Other wide-field surveys have obtained similar results on real-bogus classification by using convolutional neural networks (e.g. Duev et al. Reference Duev2019; Killestein et al. Reference Killestein2021).

3.4. Completeness test

Next, we check the completeness of Tscore classifier for real-bogus classification. For any threshold t on the score, we obtain the completeness or recovery rate as a function of magnitude, using

\[\mathrm{Completeness} = \frac{\rm TP}{\rm TP + FN},\]

where TP is the number of true positives in the test set that were correctly classified as real and FN is the number of false negatives (= positives misclassified as bogus). Figure 6 shows how the classifier performs as a function of magnitude for the two cases of asteroids and supernovae. In both cases, there is a clear trend for brighter transients to have a higher recovery rate at a given threshold. The high recovery rate of asteroids can be explained by the fact that they are over-represented in our training set compared to SN-like transients. Asteroids also appear mostly as isolated sources and are only rarely blended with galaxies, while supernovae and kilonovae are mostly blended with their host galaxy. With a threshold of $t=30$ , we attain 99.5% completeness for both asteroid (TP=55 170,170, FN=24) and supernova (TP=441, FN=2) classifications. While the SN sample has 443 detection images, these are from only 26 host galaxies and hence its statistical significance is smaller than it may seem. Moreover, there is reason to expect variations in completeness as a function of SN magnitude (at the time of detection), SN offset from the host galaxy nucleus, morphological type of the host galaxy, and redshift.

Table 3. Summary of preliminary BNS detection alerts in O3. For GW190425, we list updated information from Abbott et al. (Reference Abbott2020b)

To overcome this weakness, we test the completeness also for the DR3 SNe sample, which has fewer images but more (117 vs 26) host galaxies. Using $t=30$ again, we find high completeness ( $\sim97\%$ ; TP=148, FN=5) to the lowest magnitudes, comparable to other recently reported results (Killestein et al. Reference Killestein2021, with $\sim $ 97%). We checked the fives failures of the classifier and noticed that in three cases the FNs are superimposed on the host nucleus within $<1$ arcsec offset, so that the classification may be affected by residual features from poor subtraction. The two SNe outside the nuclei of their hosts are OGLE 16etd, a classical SN II found close to maximum light in a crowded region, and SN 2018aau, an apparently hostless SN Ia detected $\sim 2$ weeks after maximum. The latter has two observations, one classified as TP and another one as FN. The classification scores of a given transient candidate can indeed vary with seeing conditions and star density in the field, which affect image subtraction.

4.1. The search for an EM counterpart to S190425z/GW190425

On 2019 April 25 08:18:05.017 (UT), the LIGO Livingston observatory alone identified a GW chirp signal with a false alarm rate of one per 69 000 yr and a signal-to-noise ratio (SNR) of 12.9 (Abbott et al. Reference Abbott2020b). The LIGO Hanford facility was offline when the event was detected and the Virgo facility did not contribute to its detection due to a low SNR=2.5. This signal has strong evidence for the mass of one or both components to be consistent with a neutron star ( $\mathrm{HasNS}= 100\%$ ; $\mathrm{BNS}=100\%$ ). It also shows a clear time–frequency map with the characteristic upwardly sweeping chirp pattern expected for an inspiralling BNS system, which is what had been seen in GW170817 (see Figure 1 of Abbott et al. Reference Abbott2017a). The latest inferred position of this event was poorly localised in the sky, with a 90% credible sky area covering $8\,284\, \mathrm{deg}^{2}$ , but the distance inferred from the GW signal was $159^{+69}_{-71} \mathrm{Mpc}$ (see Abbott et al. Reference Abbott2020b for details). Unlike the case of GW170817, where the position constraints were $\sim 250\times$ tighter, this situation makes it difficult to find faint EM counterparts quickly. While the extremely wide-field Zwicky Transient Facility (ZTF) did search $\sim 8\,000\, \mathrm{deg}^2$ and reported two potential candidates, they needed two nights of observing to cover the area (Coughlin et al. Reference Coughlin2019). Some searches even reached sufficient depth to detect GW170817-like kilonovae, although their areal coverage was incomplete (Hosseinzadeh et al. Reference Hosseinzadeh2019).

The first optical detection by SkyMapper was made from the prioritised list of target fields that were observed about 6 h after the merger. We let the search mode (Section 2.2) run on a relatively small fraction of the southern localisation region ( $\sim125$ square degrees), acknowledging a small probability of finding the associated EM counterpart, as the area covered only $\sim1\%$ of the initial BAYESTAR map and $\sim3\%$ of its integrated probability. We found two transient candidates that had two detections separated by about 10 min, but neither of them had visible host galaxies. After the discovery of two candidates by the Zwicky Transient Facility (ZTF) survey, we moved onto a new phase to help determine whether either transient might be related to the GW signal. In Section 4.2, we perform a simple experiment to test our new software implementation with the SkyMapper observations of the BNS merger GW190425.

As part of our Target-of-Opportunity programme at the ANU 2.3m Telescope, we obtained a series of early spectra of the potential counterpart ZTF19aarykkb (Coughlin et al. Reference Coughlin2019). We used the WiFeS instrument (Dopita et al. Reference Dopita, Hart, McGregor, Oates, Bloxham and Jones2007), whose integral-field nature allows us to obtain spectra of both the transient and the host galaxy simultaneously. These could identify a kilonova from its spectral features being atypical for better known transient types (e.g. cataclysmic variable stars or supernovae) and would also provide constraints on the environment of the transient as a birth or merger place for the progenitor binary (e.g. Levan et al. Reference Levan2017). We obtained our first 850 s exposure starting at 2019 April 25 17:40:40 UT, followed by $5\times1\,200\,\mathrm{s}$ spectra, using the $R\sim3\,000$ gratings, B3000+R3000, covering the wavelength range of 320–990 nm. The WiFeS spectra showed the bright ZTF candidate (r= 18.63 mag) having spectral features similar to those of a young Type II supernova (Chang et al. Reference Chang2019b).

A number of additional EM candidates were discovered by other facilities (e.g. Hosseinzadeh et al. Reference Hosseinzadeh2019). At the location (RA=17:02:19.2, Dec=-12:29:08.2) of the Swift UVOT candidate by Breeveld et al. (Reference Breeveld2019) we did not see any transient down to an $i=20$ at the two epochs: 2019 April 25 14:17:00 and 2019 April 25 14:25:53 UT (Chang et al. Reference Chang, Wolf, Onken, Luvaul and Scott2019c). Later work suggested the Swift detection was likley to be an ultraviolet flare on an M dwarf star (Lipunov et al. Reference Lipunov2019; Bloom et al. Reference Bloom, Zucker, Schlafly, Finkbeiner, Martinez-Palomera, Goldstein and Andreoni2019). As discussed by previous studies, M dwarf flares are the most common transient phenomenon that may confuse searches for kilonova candidates from GW events (e.g. van Roestel et al. Reference van Roestel2019; Chang et al. Reference Chang, Wolf and Onken2020).

Table 4. Transient Selection for GW190425

4.2. SkyMapper search results for GW190425

To get a realistic view of the current performance of the SkyMapper transient facility, we reran the initial search mode data from the first night of the GW190425 event in O3 through the AlertSDP as described above. Each target field in the $\sim125$ square degrees of the localisation region was visited twice, separated by about 10 min. We applied the same filtering scheme for the selection of transient candidates as introduced in Section 3. Table 4 shows how each filtering step reduced the number of candidates from 41 497 to 359, which is a >100-fold reduction in the number of candidates requiring human inspection. Before the vetting process, this number was reduced even further to 9 candidates (1 in $\sim5\,000$ ) by considering only those with two detections. Each of the remaining candidates could be one of the following categories: (i) an isolated transient in a position which was not reported before, (ii) a galaxy with a SN at the time of observation, (iii) an as-yet unidentified variable object that has historical detections, or (iv) a bogus detection. We note that there was one known SN that was not detected in our search, because it was fainter than our detection limit (i = 20.4). Five of the sources were clear bogus detections caused by poor image subtraction. The remaining four candidates are persistent point sources in SMSS DR3 and seen to vary by 0.2–0.7 mag over the 5-year span of data in DR3. Thus, we conclude that none of them were related to GW190425.