Kilonova Luminosity Function Constraints Based on Zwicky Transient Facility Searches for 13 Neutron Star Merger Triggers during O3

For each of the 13 GW triggers followed up by ZTF, we systematically identified transient candidates within the localization region and ruled them out using various metrics. Below, we summarize the transient filtering process and results from our candidate vetting.

The GROWTH team has three independent database systems to retrieve interesting objects in real time: the GROWTH Marshal (Kasliwal et al. 2019a), the Kowalski ⁶² system (Duev et al. 2019), and the Alert Management, Photometry and Evaluation of Lightcurves (AMPEL) system (Soumagnac & Ofek 2018; Nordin et al. 2019). Each platform retrieves a stream of AVRO packet alerts (Patterson et al. 2019) containing significant object detections identified by the ZTF image subtraction pipeline, defined as a >5σ change in brightness relative to a reference image (Masci et al. 2018). Each of these objects undergoes a series of filtering steps in order to identify candidates that could be interesting to pursue for follow-up. The following criteria were common for all three queries.

• 1.  
  Positive subtraction. The object must have brightened relative to the reference image.
• 2.  
  Astrophysical. The object must have a real bogus (rb) score >0.25 or a deep learning (drb) score >0.8 (Duev et al. 2019; Mahabal et al. 2019) for it to be considered astrophysical.
• 3.  
  Not stellar. The object must be >2'' away from a cataloged point source in the PanSTARRS Point Source Catalog (Tachibana & Miller 2018).
• 4.  
  Far from a bright source. The object must be at least 20'' away from a bright (m_AB < 15 mag) star to avoid blooming artifacts.
• 5.  
  Not moving. The object must have at least two detections separated by at least 15 minutes to reject asteroids (moves <4'' hr⁻¹)
• 6.  
  No previous history. The object must not have any historical detections in the ZTF alert stream prior to the GW merger time.

While the GROWTH Marshal queried all fields triggered as part of the ToO search, the Kowalski and AMPEL queries searched for candidates in both serendipitous and triggered data within the 95% contour of the latest sky map that was available. The AMPEL query ⁶³ had further image quality cuts performed to reject poor subtractions based on morphology, an additional cut based on proximity to known solar system objects, and another cut based on cross-matching to the Gaia Data Release 2 (DR2) catalog and PS1 to identify likely stellar sources.

All candidates that passed the filtering criteria were saved to the GROWTH Marshal for further vetting in real time by a dedicated team of scanners. If a transient was consistent with the nucleus of a galaxy and the mid-infrared colors (based on the Wide-field Infrared Survey Explorer (WISE) catalog; Wright et al. 2010) of the host galaxy were consistent with active galactic nuclei (AGNs), the candidate was deemed unrelated.

All viable candidates were promptly announced to the worldwide community via GCN circulars, and many teams (not only GROWTH) triggered follow-up observations for many of our candidates. ⁶⁴ Using the GROWTH Marshal system, we prioritized and triggered follow-up of candidates that exhibited rapid photometric evolution (faster than 0.3 mag day⁻¹), showed red colors, or were close to a host galaxy with a redshift consistent with the GW distance constraint.

We now briefly describe how we ruled out the association between vetted counterpart candidates and the GW event. A detailed account of every candidate announced via GCN is in Appendix B.

The GROWTH team obtained follow-up with the following facilities to characterize the photometric and/or spectroscopic evolution: the Liverpool Telescope (LT; Steele et al. 2004), Lowell Discovery Telescope (LDT, ⁶⁵ formerly known as the Discovery Channel Telescope), Las Cumbres Observatory (LCO; Brown et al. 2013), Apache Point Observatory (APO; Huehnerhoff et al. 2016), Kitt Peak EMCCD Demonstrator (KPED; Coughlin et al. 2019b), Lulin One-meter Telescope (LOT; Huang et al. 2005), GROWTH-India telescope (GIT; ⁶⁶ V. Bhalerao et al. 2020, in preparation), Palomar 60 inch telescope (P60; Cenko et al. 2006), Palomar 200 inch Hale Telescope ⁶⁷ (P200), Keck Observatory, ⁶⁸ Gemini Observatory, ⁶⁹ Southern African Large Telescope ⁷⁰ (SALT), Himalayan Chandra Telescope ⁷¹ (HCT), and Gran Telescopio Canarias ⁷² (GTC). Figures 6 and 7 illustrate examples of follow-up by the GROWTH team on some ZTF candidates. The specific instrument configurations and data reduction methods are described in Appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The follow-up observations include both photometric and spectroscopic data. Moreover, the association of a candidate with a GW trigger was rejected if its properties fell into one or more of the categories described as follows.

• 1.  
  Inconsistent spectroscopic classification. We ruled out candidates that could be spectroscopically classified as SNe, AGNs, cataclysmic variables (CVs), and other flare stars. We used SNID (Blondin & Tonry 2007) and dash (Muthukrishna et al. 2019) to classify the SNe and AGNs found in our searches. The CVs and variable stars often showed hydrogen features at zero redshift.
• 2.  
  Inconsistent distance. We ruled out candidates whose spectroscopic redshift was not consistent with the GW distance within 2σ. We cross-matched the transient positions with the Census of the Local Universe (CLU; Cook et al. 2019) galaxy catalog and the NASA Extragalactic Database (NED) to look up host redshifts where available. We also cross-matched the candidates against the Photometric Redshifts Legacy Survey (PRLS; Zhou et al. 2020) catalog and reported the photometric redshifts when the spectroscopic redshift was unavailable.
• 3.  
  Slow photometric evolution. As kilonovae are expected to evolve faster than SNe, we ruled out candidates that evolved slower than 0.3 mag day⁻¹. We used ForcePhot ⁷³ (Yao et al. 2019), a forced photometry package, to examine the transient light curves. To quantify the evolution of a given transient, we define the parameter α_f = Δm/Δt [mag day⁻¹], where f corresponds to the filter used to determine the variation in magnitude (Δm) over time (Δt). A positive α indicates a fading source, while a negative α describes a rising source. The baseline (Δt) is defined to be the number of days it takes an object to rise from its discovery to its peak magnitude (α < 0) or the number of days it takes the transient to fade from peak to undetectable by ZTF (α > 0). We used a minimum time baseline of 3 days to compute slopes.
• 4.  
  Outside of the latest LALInference map. The majority of the candidates were selected and announced via GCN based on the promptly available BAYESTAR map (Singer & Price 2016). When the LALInference map was made available, if a candidate was outside the 90% probability contour, we rejected it.
• 5.  
  Artifacts. Most of the ZTF ghosts and artifacts are well known (Bellm et al. 2018; Masci et al. 2018) ⁷⁴ and masked automatically. Additionally, we take further precautions by ignoring transients close to bright stars in our initial vetting. However, for example, our extensive analysis revealed a subtle gain mismatch in the reference images that posed as a faint and fast transient (see discussion related to ZTF19aassfws in Appendix B). All references for ToOs were rebuilt after this artifact was identified.
• 6.  
  Asteroids. Sometimes slow-moving asteroids, especially near stationary points, can mimic a fast-fading transient (Jedicke et al. 2016). For these objects, either a more careful inspection of the centroids or movement in follow-up imaging served as the reason for rejection.
• 7.  
  Previous activity. Candidates were rejected if they showed previous detections prior to the GW merger time in other surveys, e.g., the Catalina Real Time Survey (CRTS; Djorgovski et al. 2011), Palomar Transient Factory (PTF; Law et al. 2009), intermediate Palomar Transient Factory (iPTF; Cao et al. 2016; Masci et al. 2017), and PS1 (Tachibana & Miller 2018).

Some candidates prompted panchromatic follow-up. We followed up five candidates in the ultraviolet and X-ray with the Neil Gehrels Swift Observatory (see the Appendix for details). We followed up two candidates in the radio with the Arcminute Microkelvin Imager (AMI) and one with the Karl G. Jansky Very Large Array (VLA; see the Appendix for details). All candidates, grouped by GW trigger, are listed in Tables 2–10, along with their respective rejection criteria.

We complemented our real-time analysis described above with a deeper offline search by relaxing the selection criteria (e.g., requiring only one detection instead of two). The following steps describe our offline search.

• 1.  
  We used Kowalski to query the ZTF database looking for any source (i) located within 95% of the most updated sky map, (ii) never detected before the merger time, (iii) with at least one detection within 72 hr of merger, (iv) with the last detection occurring within 10 days of the first detection, and (v) passing real/bogus thresholds of rb > 0.5 and drb > 0.8 (or braai > 0.8; Duev et al. 2019). Further details on the selection criteria will be described in I. Andreoni et al. (2020, in preparation).
• 2.  
  Forced point-spread function (PSF) photometry was performed at the location of each transient candidate using ForcePhot, setting a detection threshold of signal-to-noise ratio (S/N) > 3, where the images were available.
• 3.  
  The flux measured using forced photometry was stacked nightly in each band, allowing us to become sensitive to fainter sources when multiple images were available on the same night.
• 4.  
  The rising and fading rates were computed in each band with a linear fit before and after the brightest data point of each light curve. A time baseline of >3 hr was required for the fit to be performed.
• 5.  
  Candidates were selected with a fading rate more rapid than 0.3 mag day⁻¹ or rising rate faster than 1 mag day⁻¹. We rejected candidates still detected after 6, 12, and 14 days after the merger time in the g, r, and i bands, respectively. More details are given in I. Andreoni et al. (2020, in preparation).

The Kowalski query initially returned 8026 sources for the 13 GW triggers. Applying all of the selection criteria described above, 453 candidates survived the automatic cuts. Of these, 21 had at least two ZTF alerts, and 432 had only one ZTF alert (additional detections were recovered by forced photometry and stacking).

Of the 21 sources with at least two detections in the ZTF alert stream, only five candidates passed visual inspection of the images and light curves: ZTF19acbxacj was an AGN candidate (Assef et al. 2018; Bailer-Jones et al. 2019); ZTF19abwsfsl was a cataloged CV (Gaia Collaboration 2018); ZTF19acbqtue was followed up with the Gemini Multi-Object Spectrograph (GMOS-N), and a quiescent source was found at g = 24.69 ± 0.07 mag with a color g − i = 1.89 mag, consistent with an M dwarf (West et al. 2011); ZTF19abyndjf was a fast-evolving transient without an obvious host galaxy; and ZTF19acbwtmt was hostless and had a previous detection in the PS1-DR2 catalog from 2012 (see Figure 8). For the last two candidates, upper limits between the GW merger time and the first transient detection disfavor their multimessenger association with S190930t or S190910d, respectively (see Figure 8).

Zoom In Zoom Out Reset image size

Of the 432 sources with only one detection in the ZTF alert stream, only nine candidates passed the visual inspection of the images and light curves. Most other candidates were ruled out as stellar flares, image subtraction artifacts, asteroids, or sporadic nuclear variability. Of these nine candidates, six had photometric or spectroscopic redshifts of the host galaxy too far to be consistent with the GW distance. All three remaining candidates were found during follow-up of S190901ap: ZTF19abvpeir, ZTF19abvozxv, and ZTF19abvphxm. All three are likely SNe or AGNs, given that their absolute magnitudes at the distance of their putative hosts are between −18.0 and −18.8 mag and their locations are consistent with the Galaxy nuclei. We show some light curves and host galaxies in Figure 8.

In summary, all candidates were ruled out as possible kilonovae in both the real-time and offline analyses.

We estimate the joint probability of zero kilonova detections as a product of (1 − q_i ) terms, where q_i is the enclosed probability for each event as listed in Table 1 (see the sixth column; we used LALInference probability where available). If we were sufficiently sensitive to finding kilonovae in all 13 GW events, the joint probability of zero detections would be only 0.017%. However, each merger had a different observed depth, observed cadence, and GW distance estimate and thus a different sensitivity to detecting kilonovae.

First, we use the median image depth for each trigger and the median GW distance to each trigger to compute a median absolute magnitude sensitivity limit. We correct the median absolute magnitude for the median extinction along the line of sight. In Figure 9, in each luminosity bin, we compute the joint nondetection probability only for the subset of events for which the ZTF observations were sufficiently sensitive. We find that ZTF follow-up of four (six) GW events had a sensitivity deeper than −16.0 mag (−16.6 mag), and the joint nondetection probability is only 4.5% (0.34%). Moreover, three of the four (four of the six) had a preliminary BNS classification, and for all three, the ZTF follow-up began within 4 hr of merger (see Table 1).

Zoom In Zoom Out Reset image size

Second, we use injection and recovery of fake sources to better quantify both the degree of variation in the depths of individual exposures and the spatial variation in the GW distance estimates. We use an open-source tool called simsurvey ⁷⁵ (Feindt et al. 2019). As input, the tool takes a list of ZTF pointings (observation time, R.A., decl., limiting magnitude, filters, processing success of each CCD quadrant). We inject 10,000 sources distributed according to the 3D GW sky map probability distribution in each luminosity bin (50 bins between −10 and −20 mag). Initially, we assume that each kilonova stays at a constant luminosity between the merger time and 3 days after merger. We require a single observation at the necessary depth for recovery. In addition to losing sources in unobserved fields, we lose sources that land in ZTF chip gaps, chips that failed processing, or chips that were less sensitive due to higher line-of-sight Galactic extinction. This tool does not take into account any detections that would be lost to inefficiency in the software pipeline.

The recovery fraction for each event is shown in Figure 10. We convert this to a joint nondetection probability estimate by multiplying (1 − _pi ) in each luminosity bin and overlay this as discrete points on the median estimates above in Figure 9. We find consistent results; the joint probability of zero detections at −16.1 mag (−16.6 mag) is only 4.2% (0.8%). If we separate the NSBH and BNS populations, the joint nondetection probability at −16.6 mag is 9.7% for NSBH and 7.9% for BNS. This is not surprising, as the BNS triggers were, on average, closer than the NSBH triggers. We note that this application of simsurvey is different compared to previous applications for SN rates, which were uniform in volume. Taking into account the exact 3D GW sky map is more accurately representative of our success in searching for the counterpart to a GW source on an event-by-event basis.

Zoom In Zoom Out Reset image size

Third, in addition to spatial variations in depth and distance, we take into account the possible time variations in the light curves of kilonovae (Figure 11). The time window for our observations is limited to within 3 days of merger. We relax the constant luminosity assumption above and inject kilonovae into simsurvey that fade linearly between zero and 1 mag day⁻¹. In Figure 11, we color-code the recovery efficiency for a given peak luminosity and photometric decay rate in any filter (g or r band). Any slice of this plot can be converted to a joint probability of zero detections as a function of absolute magnitude. We compare to the GW170817 benchmark of an extrapolated peak of −16.6 mag day⁻¹ and a fade rate of 1 mag day⁻¹. We find a joint probability of zero detections of 7%.

Zoom In Zoom Out Reset image size

Fourth, in addition to spatial and time variations, we inject kilonova models into simsurvey and calculate the recovery fraction. We use the best fit to GW170817 from the kilonova model grid in Dietrich et al. (2020) computed using the radiative transfer code POSSIS (Bulla 2019). This model fit assumes a rise time, color evolution, and viewing angle of GW170817. The joint nondetection probability is 30%. Even if all kilonova ejecta parameters were similar to GW170817, the viewing angle would be different for different events. Assuming random viewing angles drawn from a distribution uniform in cos(θ), we inject a model grid and find that the joint nondetection probability is 49%. We caution that the model used here underestimates the early g-band flux of GW170817 by 0.3 mag; thus, the recovery fraction estimated here could also be underestimated.

Next, we consider the implications of the zero-detection probability function on the underlying luminosity function. Let us say the luminosity function dictates that a fraction f_b of kilonovae are brighter than a given absolute magnitude. Then,

$$\begin{eqnarray*}&&(1-\mathrm{CL})=\displaystyle \prod _{i=1}^{N}(1-{f}_{b}^{* }{p}_{i}),\end{eqnarray*}$$

where CL is the confidence level and _pi is the event-by-event probability of detection. At a given absolute magnitude, we compute _pi as the recovery fraction from the simsurvey injections for the fading and flat light-curve evolution estimates discussed above that take into account the spatial variation in distance, depth, and enclosed probability. Solving for f_b at a 90% confidence level, we plot our results in Figure 12. At the bright end, we find that no more than ≈40% of kilonovae can be brighter than −18 mag. At the faint end, our observations place no constraints on the luminosity function fainter than −15.5 mag. The luminosity of GW170817 at the merger time is unknown, and various models predict diverse rates of evolution in that first day after merger. As discussed above, we use an extrapolated peak of −16.6 mag and a fade rate of 1 mag day⁻¹ for GW170817 as a benchmark. We find that no more than <57% (<89%) of kilonovae could be brighter than −16.6 mag for the flat (fading) light-curve assumptions.

Zoom In Zoom Out Reset image size

The GW triggers had a very wide range of FARs. Weighting by the available low-latency values for the terrestrial probability (t_i ), we fold this into our luminosity function constraint as

$$\begin{eqnarray*}&&(1-\mathrm{CL})=\displaystyle \prod _{i=1}^{N}(1-{f}_{b}^{* }{p}_{i}^{* }(1-{t}_{i})).\end{eqnarray*}$$

In Figure 12, we show that the resulting constraints on f_b (red line) are worse only by a difference of ≈10%.

Next, we investigate the implications of this constraint on the kilonova parameter space. There are no theoretical luminosity functions available in the literature that we can directly compare to. A model grid is available (Kasen et al. 2017) as a function of three parameters: ejecta mass M_ej, ejecta velocity v_ej, and lanthanide fraction X_lan. The best-fit model to GW170817 from Kasen et al. (2017) suggested two components: a blue kilonova (0.025 M_⊙, 0.3 c, 10^−4.5) and a red kilonova (0.04 M_⊙, 0.1 c, 10⁻²). The blue component dominates at an early time and is more relevant to the ZTF searches described in this paper. Comparing to our luminosity function constraints, we find that our limits suggest that a wide range of parameters are allowed, e.g., M_ej = [0.03, 0.1] M_⊙, v_ej = [0.05, 0.3] c, and X_lan = [10⁻⁵, 10⁻⁴]; stricter distributions that yield a brighter kilonova population (e.g., higher ejecta mass or lower lanthanide fraction) are not allowed. Thus, some kilonovae must have M_ej ≤ 0.03 M_⊙ or X_lan > 10⁻⁴ to be consistent with the ZTF constraints.

Similarly, we compare our luminosity function constraints to the kilonova grid from Dietrich et al. (2020) computed using the radiative transfer code POSSIS (Bulla 2019). In addition to the observer viewing angle, this grid depends on three parameters: the dynamical ejecta mass (M_ej,dyn), the postmerger wind ejecta mass (M_ej,pm), and the half-opening angle of the lanthanide-rich ejecta component (ϕ). A model with M_ej,dyn = 0.005 M_⊙, M_ej,pm = 0.05 M_⊙, and ϕ = 30° provides a good fit to GW170817 (see Figure 8 of Dietrich et al. 2020). As shown in Figure 12, our constraints suggest that some kilonovae must be fainter than GW170817, i.e., must have either M_ej,dyn < 0.005 M_⊙, M_ej,pm < 0.05 M_⊙, or ϕ > 30°.

We used the 1 and 2 m telescopes available at the LCO global network to follow up sources discovered with the ZTF. The images were taken with the Sinistro and Spectral cameras (Brown et al. 2013) at the 1 and 2 m, respectively, and scheduled through the LCO Observation Portal. ⁷⁶ The exposure time varied depending on the brightness of the object, yet our requests would normally involve three sets of 300 s in the g and r bands. After stacking the reduced images, we extracted sources using the SExtractor package (Bertin & Arnouts 2010), and we calibrated magnitudes against PS1 (Chambers et al. 2016) objects in the vicinity. For nuclear transients located <8'' from their potential host, we use the High Order Transform of Psf ANd Template Subtraction code (HOTPANTS; Becker 2015) to subtract a PSF-scaled PS1 template previously aligned using SCAMP (Bertin 2006). The photometry for the nuclear candidates follows the same procedure described before but in the residual image. The images obtained with the LT were acquired using the IO:O camera with the Sloan griz filter set. They were reduced using the automated pipeline, which performs bias subtraction, trimming of the overscan regions, and flat-fielding. The image subtraction takes place once a PS1 template is aligned, and the final data come from the analysis of the subtracted image.

We used the Electronic Multiplier CCD camera at KPED to take hour-long exposures in the r band to follow up candidates. After stacking the images and following standard reduction techniques, we calibrate the extracted sources using PS1 sources in the field. When the candidate has a host galaxy, we perform image subtraction as described for the LCO.

We obtained data with the GMOS-N (Allington-Smith et al. 2002; Hook et al. 2004; Gimeno et al. 2016), mounted on the Gemini-North 8 m telescope on Maunakea. Data were analyzed after stacking four 200 s exposures in the g and i bands. The reductions were performed using the python package DRAGONS ⁷⁷ provided by the Gemini Observatory. We used PS1 sources in the field to calibrate the data.

We used the LOT at the Lulin Observatory in Taiwan to follow up candidates discovered with the ZTF. The standard observations involved 240 s in the g', r', and i' bands. The reduction followed standard methods, and the sources were calibrated against the PS1 catalog. No further image subtraction was applied to the images acquired with the LOT.

We used the 0.7 m robotic GIT equipped with a 4096 × 4108 pixel back-illuminated Andor camera for LVC event follow-up during O3. The GIT is situated at the IAO (Hanle, Ladakh). We used both tiled and targeted modes for the follow-up for different GW triggers. Tiled observations typically comprise a series of 600 s exposures in the Sloan Digital Sky Survey (SDSS) r' filter. Targeted observations were conducted with varying exposure times in the SDSS u', g', r', and i' filters. All data were downloaded in real time and processed with the automated GIT pipeline. Zero-points for photometry were calculated using the PanSTARRS catalog (Flewelling 2018), downloaded from Vizier. The PSF photometry was performed with PSFEx (Bertin 2011). For sources with significant host background, we performed image subtraction with pyzogy (Guevel & Hosseinzadeh 2017), based on the ZOGY algorithm (Zackay et al. 2016).

Additionally, we obtained photometric data with the Spectral Energy Distribution Machine (SEDM; Blagorodnova et al. 2018; Rigault et al. 2019) on the P60 telescope. The processing is automated and can be triggered from the GROWTH Marshal. Standard requests involved g-, r-, and i-band imaging with the Rainbow Camera on the SEDM in 300 s exposures. The data are later reduced using a python-based pipeline that applies standard reduction techniques and a customized version of the Fremling Automated Pipeline (FPipe; Fremling et al. 2016) for image subtraction.

We used the imaging capabilities of the OSIRIS (Cepa et al. 2005) camera at the GTC to obtain 60 s exposures in the r band. Standard reduction techniques were applied to the data, and we used PS1 sources to calibrate the flux.

We obtained follow-up imaging of candidates with the Wafer Scale Imager for Prime (WASP) and the Wide-field Infrared Camera (WIRC; Wilson et al. 2003), both on the P200 telescope. For WASP data, a python-based pipeline applied standard optical reduction techniques (as described in De et al. 2020a), and the photometric calibration was obtained against PS1 sources in the field. The WIRC data were treated similarly using the same pipeline, but they were additionally stacked using Swarp (Bertin et al. 2002), while the calibration was done using the Two Micron All Sky Survey point-source catalog (Skrutskie et al. 2006).

We obtained imaging of one candidate using the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995), mounted at the Keck I telescope. Our data were taken in the g and i bands reaching m_AB ≈ 24. The data were reduced following standard methods.

We used the Large Monolithic Imager (LMI; Massey et al. 2013) on the 4.3 m LDT at Happy Jack, Arizona, to follow up the ZTF discoveries. Observations were conducted with the SDSS r filter for 90 s each, and the data were reduced using the photopipe ⁷⁸ pipeline. The magnitudes were calibrated against the SDSS or Gaia catalogs (Ahumada et al. 2020) using the conversion scheme provided in Gaia documentation. ⁷⁹

We used the Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005) mounted on the Neil Gehrels Swift Observatory (hereafter referred to as Swift; Gehrels et al. 2004) to follow up interesting sources and track down their UV evolution. The ToO observations were scheduled in the v, b, u, w1, m2, and w2 bands for an average of 320 s exposure^–1. We used the products of the Swift pipeline to determine the magnitudes. ⁸⁰

We observed candidate counterparts of S200213t using the Astrophysical Research Consortium Telescope Imaging Camera (ARCTIC; Huehnerhoff et al. 2016) at the APO 3.5 m. We obtained dithered 120 s exposures binned 2 × 2 in the u, g, r, i, and z bands. Images were bias-corrected, flat-fielded, and combined using standard IRAF packages (noao, imred, and ccdred). SExtractor (Bertin & Arnouts 2010) was used to find and photometer point sources in the images using PSF photometry, and a photometric calibration to PanSTARRS field stars was performed (without filter corrections).

All photometry presented in the light curves and tables in this paper is corrected for galactic extinction using dust maps from Schlafly & Finkbeiner (2011).

We observed the field of ZTF19aassfws with the VLA in its B configuration on 2019 May 10 starting at 07:19:15 UT and 2019 June 4 starting at 08:20:32 UT. Our observations were carried out at a nominal central frequency of 3 GHz. We used 3C 286 as our bandpass and absolute flux calibrator and J1927+6117 as our complex gain calibrator. Data were calibrated using the standard VLA automated calibration pipeline available in the Common Astronomy Software Applications (CASA) package. We then inspected the data for further flagging and imaged interactively using the CLEAN algorithm. The image rms was ≈5.2 μJy for the first epoch and ≈4.6 μJy for the second epoch. Within a circular region centered on the optical position of ZTF19aassfws and of radius ≈21 (comparable to the nominal half-power beamwidth of the VLA at 3 GHz and for the B configuration), we find no significant radio emission. Thus, we set upper limits on the corresponding 3 GHz flux density of ≲16 and ≲14 μJy, respectively, for the first and second epochs.

Using the GROWTH Marshal, we regularly triggered the Liverpool Telescope Spectrograph for the Rapid Acquisition of Transients (SPRAT; Piascik et al. 2014). SPRAT uses a 18 slit, which provides a resolution of R = 350 at the center of the spectrum. The data were reduced using the automated pipeline, which removes low-level instrumental signatures and then performs source extraction, sky subtraction, and wavelength and flux calibration.

We observed a number of transient candidates during classical observing runs with the P200 Double Spectrograph (DBSP) during O3. For the setup configuration, we used 10, 15, and 2'' slit masks; a D55 dichroic; a blue grating of 600/4000; and a red grating of 316/7500. Using a custom PyRAF DBSP reduction pipeline (Bellm & Sesar 2016), ⁸¹ we reduced our data.

We obtained several optical spectra with the 10.4 m GTC telescope (equipped with OSIRIS). We used the R1000B and R500R grisms for our observations, typically using a slit of width of 12. We used standard routines from the Image Reduction and Analysis Facility (IRAF) to perform our data reduction.

We observed ZTF19aarykkb using the DeVeny spectrograph mounted on the 4.3 m LDT. We obtained 22.5 minute exposures at an average airmass of 1.5. We used the DV2 grating (300 g mm⁻¹, 4000 Å blaze) for this observation. Our spectra cover a wavelength range of approximately 3600–8000 Å.

In addition, we obtained a spectrum of ZTF20aarzaod with SALT (Buckley et al. 2003), using the Robert Stobie Spectrograph (RSS; Burgh et al. 2003), covering a wavelength range of 470–760 nm with a spectral resolution of R = 400. We triggered special GW follow-up program 2018-2-GWE-002 and reduced the data with a custom pipeline based on PyRAF routines and the PySALT package (Crawford et al. 2010).

Low-resolution spectra using the 2 m HCT were obtained using the HFOSC instrument; ZTF19aarykkb was observed using grisms Gr7 (3500–7800 Å) and Gr8 (5200–9000 Å), while AT2019wxt was observed using Gr7. The spectra were bias-subtracted, cosmic rays were removed, and the 1D spectra were extracted using the optimal extraction method. Wavelength calibration was effected using the arc lamp spectra FeAr (Gr7) and FeNe (Gr8). Instrumental response curves generated using spectrophotometric standards observed during the same night were used to calibrate the spectra onto a relative flux scale. The flux-calibrated spectra of ZTF19aarykkb from the two grisms were combined to a single spectrum covering the wavelength range 4000–9000 Å.

We obtained spectroscopy with the GMOS-N, mounted on the Gemini-North 8 m telescope on Maunakea, by combining six 450 s exposures on the R400 and B600 gratings. We used the GMOS long-slit capability and reduced the data following standard PyRAF techniques.

We obtained near-infrared spectroscopy of candidates using NIRES on the Keck II telescope. The data were acquired using standard ABBA dither patterns on the target source, followed by observations of an A0 telluric standard star close to the science target. The spectral traces were extracted using the spextool package (Cushing et al. 2004) for both the science target and standard star. The final spectra presented here were stacked from all individual dithers, followed by flux calibration and telluric correction using the xtellcor package (Vacca et al. 2003).

We obtained spectra using the LRIS on the Keck I telescope. The 600/4000 grism was used on the blue side and the 600/7500 grating on the red side, providing wavelength coverage between 3139–5642 Å (blue) and 6236–9516 Å (red). The exposure time was 600 s on both sides. The spectrum was reduced using LPipe (Perley 2019) with BD+28 as a flux calibrator. The red and blue relative fluxes are scaled by matching synthetic photometry to colors inferred from photometry of the transient.

For candidates identified within the sky map of GW190425, see Coughlin et al. (2019a). Two candidate counterparts of GW190425z, ZTF19aarykkb and ZTF19aarzoad, were observed with the AMI Large Array at 15 GHz on 2019 April 26 (Rhodes et al. 2019). No radio emission was found to be associated with any of these candidates.

We summarize the candidate counterparts to S190426c in Table 2 and the follow-up photometry in Table 11. Next, we discuss why we conclude that each one is unrelated.

B.2.1. Spectroscopically Classified

B.2.2. Slow Photometric Evolution

B.2.3. Stellar

B.2.4. Artifacts

No candidates were identified in the ZTF follow-up of the small localization of GW190814.

We summarize the candidate counterparts to S190901ap in Table 3 and the follow-up photometry in Table 12. Next, we discuss why we conclude that each one is unrelated.

B.4.1. Spectroscopically Classified

B.4.2. Slow Photometric Evolution

We summarize the candidate counterparts to S190910d in Table 4 and the follow-up photometry in Table 13. Next, we discuss why we conclude that each one is unrelated.

B.5.1. Spectroscopically Classified

We summarize the candidate counterparts to S190910h in Table 5 and the follow-up photometry in Table 14. Next, we discuss why we conclude that each one is unrelated.

B.6.1. Spectroscopically Classified

B.6.2. Slow Photometric Evolution

B.6.3. Artifacts

We summarize one candidate counterpart to S190923y in Table 6. Despite the small sky localization, the position of S190923y on the sky made it particularly challenging to access. For that reason, we chose to conduct a fully serendipitous search in ZTF data.

ZTF19acbmopl/AT2019rob—We found this transient with a photometric redshift of ≲0.03, consistent with the LVC distance reported, slightly off the nucleus of its host galaxy. It showed a slow evolution in both the r and g bands: α_r  = 0.03 and α_g  = 0.03.

We summarize the candidate counterparts to S190930t in Table 7 and the follow-up photometry in Table 15. Next, we discuss why we conclude that each one is unrelated.

B.8.1. Spectroscopically Classified

We summarize the candidate counterparts to S191205ah in Table 8 and follow-up photometry in Table 16. Next, we discuss why we conclude that each one is unrelated.

B.9.1. Spectroscopically Classified

B.9.2. Slow Photometric Evolution

We summarize the candidate counterparts to S191213g in Table 9 and the follow-up photometry in Table 17. Next, we discuss why we conclude that each one is unrelated.

B.10.1. Spectroscopically Classified

B.10.2. Slow Photometric Evolution

B.10.3. Artifacts

B.10.4. Stellar Sources

For candidates identified within the sky map of S200105ae and S200115j, see Anand et al. (2020).

We summarize the candidate counterparts to S200213t in Table 10 and the follow-up photometry in Table 18. Next, we discuss why we conclude that each one is unrelated. All of the transients described for this event (S200213t) were reported in Kasliwal et al. (2020).

B.12.1. Spectroscopically Classified

B.12.2. Stellar

B.12.3. Slow Photometric Evolution

B.12.4. Outside the GW Map