Constraining the Kilonova Rate with Zwicky Transient Facility Searches Independent of Gravitational Wave and Short Gamma-Ray Burst Triggers

The study of the dynamic sky took a twist when multi-messenger discoveries became reality with the discovery of astronomical transients associated with neutrinos (Hirata et al. 1987; IceCube Collaboration et al. 2018; Stein et al. 2020) or gravitational waves (GWs; Abbott et al. 2017d). The identification of an electromagnetic (EM) counterpart to the GW event GW170817 (Abbott et al. 2017d) spectacularly confirmed the predictions that short gamma-ray bursts (GRBs; e.g., Blinnikov et al. 1984; Paczynski 1986; Eichler et al. 1989; Narayan et al. 1992; Fong & Berger 2013) and fast optical/near-infrared transients called "kilonovae" or "macronovae" (KNe; e.g., Li & Paczyński 1998; Kulkarni 2005; Rosswog 2005; Metzger et al. 2010; Metzger & Berger 2012; Barnes & Kasen 2013; Tanaka & Hotokezaka 2013) can be generated by binary neutron star (BNS) or neutron star–black hole (NSBH) mergers (see also Foucart 2012; Hotokezaka et al. 2013). Launched by neutron star disruption, neutron-rich ejecta host rapid neutron capture (r-process) nucleosynthesis, leading to the production of heavy elements. The radioactive decay of such unstable heavy nuclei powers the optical/infrared KN. Recent reviews on EM emission from neutron star mergers can be found, for example, in Metzger (2019), Nakar (2019).

The observation of optical/infrared KNe is particularly valuable when it happens concurrently with the discovery of a compact binary merger in GWs. The nearby BNS merger GW170817 was confidently associated with the optical transient AT2017gfo (Coulter et al. 2017) and the multi-messenger discovery has led to hundreds of studies addressing, for example, astrophysics of energetic phenomena (e.g., Kasliwal et al. 2017; D'Avanzo et al. 2018; Margutti et al. 2018; Beniamini et al. 2019; Salafia & Giacomazzo 2020), fundamental physics (e.g., Abbott et al. 2019; Coughlin et al. 2019b), and cosmology (e.g., Abbott et al. 2017a; Hotokezaka et al. 2019). It is relevant for this work to note that even the optical/infrared photometric data alone carried evidence of heavy-element nucleosynthesis in the GW170817 merger ejecta (Andreoni et al. 2017; Arcavi et al. 2017; Chornock et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Pian et al. 2017; Shappee et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017).

The discovery of more KNe in optical surveys would allow us to estimate the production of heavy elements during neutron star mergers and better constrain the neutron star merger rate (see Section 5), therefore providing tools to tackle open questions such as whether neutron star mergers are "the site" or just "a site" for heavy-element production in the universe. This topic is discussed, for example, by Rosswog et al. (2018), Kasliwal et al. (2019), Siegel et al. (2019). The possibility of using KNe as standardizable candles for cosmology (without any GW information) is also an intriguing perspective (Coughlin et al. 2020) that can be tested when a significant population of KNe is unveiled.

Observationally, KNe can show an early blue component likely brighter toward polar angles for the first 2–3 days after merger, followed by a redder component generated from the tidal dynamical ejecta and post-merger ejecta (e.g., Cowperthwaite et al. 2017; Kasen et al. 2017; Kawaguchi et al. 2018; Bulla et al. 2019). Regardless of the viewing angle, KNe are expected to appear dimmer and fade more quickly relative to most supernovae at optical wavelengths. Supernovae typically peak at absolute magnitudes between $\sim -17$ and $\sim -19$ and, in most cases, they fade at a rate slower than 0.3 mag day⁻¹. AT2017gfo faded at a rate $\gt 0.5$ mag day⁻¹ in optical bands and it was measured to have an absolute magnitude of ${M}_{r}\sim -16$ shortly after peak. Such characteristics make it more difficult to discover KNe relative to supernovae, which are brighter at peak and remain luminous for months. Moreover, recent GW observations indicate a rate for BNS mergers of 250–2810 Gpc⁻³ yr⁻¹ (Abbott et al. 2020), which suggests that at most a few events per year are expected to occur at distances close enough for KNe to be detectable in 1 m class optical surveys. In comparison, the rate of unobscured supernovae in the nearby universe is $\sim 28.7\,\times {10}^{3}$ Gpc⁻³ yr⁻¹ (Smartt et al. 2009).

Several KN candidates were found during the follow-up of short GRBs (Perley et al. 2009; Berger et al. 2013; Tanvir et al. 2013; Gao et al. 2015; Jin et al. 2015, 2016, 2020; Yang et al. 2015; Troja et al. 2018; Lamb et al. 2019) as optical/infrared excesses on top of the GRB afterglow and one was found during the follow-up of GW170817/GRB170817A (e.g., Coulter et al. 2017). Although a marginal candidate was recently identified (McBrien et al. 2020), no confirmed KN has been found during optical surveys independently of GW triggers, to date.

Searches in optical time-domain survey data sets are complementary to GW follow-up searches. The third Advanced LIGO/Virgo observing run (O3) took place entirely within the time frame in which data used in this work were taken and it was suspended without yielding any optical (or multi-wavelength) counterpart (Andreoni et al. 2019b, 2020b; Coughlin et al. 2019a; Goldstein et al. 2019; Gomez et al. 2019; Hosseinzadeh et al. 2019; Lundquist et al. 2019; Ackley et al. 2020; Antier et al. 2020; Gompertz et al. 2020; Vieira et al. 2020; Watson et al. 2020). In addition, it is possible that neutron star mergers occurred in the nearby universe during O3 and were missed by the GW detectors because of instrument downtime and their unisotropic antenna pattern.

Motivated by this, we use the Zwicky Transient Facility (ZTF; Bellm et al. 2019b; Graham et al. 2019) at Palomar Observatory, which surveys the northern sky primarily in g and r band (i band being occasionally used for specific projects) with an instantaneous field of view of 47 deg², a pixel scale of $1^{\prime\prime}$ pixel⁻¹, and a typical limiting magnitude of $r\sim 20.5$ in 30 s exposure time, to search for serendipitous KNe. This work is based on 239,754 images acquired during the first 23 months of ZTF operations.

The paper is organized as follows: methods used to search the ZTF alert database, to perform photometry, and to select KN candidates are described in Section 2; the results of our searches are presented in Section 3; implications of our results for KN rates are discussed in Section 5. Section 6 completes the manuscript with a short summary and concluding remarks. A Planck Collaboration et al. (2016) cosmology is used throughout the paper.

In this section, we describe the data set that we searched for KNe, the data mining techniques that we adopted, and our novel scheme to vet candidates.

2.1. Data Set

ZTF conducts a public Survey of the Northern Sky for 40% of the time allocation. Public data are acquired as g, r filter pairs with a preferred cadence of three nights, with some fields in the Galactic plane or in Transiting Exoplanet Survey Satellite (Ricker et al. 2015) fields imaged with nightly cadence (for more information on ZTF goals and planning see Bellm et al. (2019a). The remaining time allocation is divided between ZTF Partnership and Caltech surveys, usually observing a smaller sky area (including the Galactic plane) at higher cadence. Data from all programs were used in our analysis. Images are processed with the ZTF real-time reduction and image-subtraction pipeline at the Infrared Processing & Analysis Center (Masci et al. 2019), using the ZOGY algorithm for image subtraction (Zackay et al. 2016). Upon source detection, alert packets are generated and serialized in Apache Avro format, and finally distributed to community brokers (Patterson et al. 2019).

This work aims to discover extragalactic transients, so we limit the search for KNe to 514 fields (47 deg² each) with low Galactic extinction by requiring $E(B-V)\lt 0.3$ mag at the central coordinates for the field based on Planck Collaboration et al. (2014) dust maps. The central coordinates of ZTF fields belong to two fixed grids, where the secondary grid is shifted from the primary grid in order to cover the chip gaps. We only take fields in the well-sampled ZTF primary grid, given that the secondary grid has significantly fewer observations. The effective area that we probed extended across ∼24,150 deg² (Figure 1). The resulting data set utilizes 239,754 images, 89,961 of which were acquired in g band, 137,403 in r band, and 12,390 in i band.

Zoom In Zoom Out Reset image size

2.2. Alert Database Query

2.3. Photometry and Candidate Selection Criteria

Forced point-spread function (PSF) photometry was performed at the location of each candidate using ForcePhotZTF (Yao et al. 2019). The median coordinates recorded in the available alerts were used to improve the location accuracy. A signal-to-noise ratio of S/N = 3 served as threshold for the photometric detection of a source.

First, we ran forced photometry on images taken starting a week before the first detection and until 2 weeks after the last detection of each candidate, in order to be less biased by spurious detections consistent with Gaussian noise in several thousands of images by considering the whole survey. Then, we stacked the forced photometry data points nightly in flux space to increase our sensitivity to faint sources. Finally, we ran forced photometry at the location of selected candidates for the whole duration of the survey and we rejected those candidates showing repeated activity that was not significant enough to be recorded in ZTF alerts. An example of how forced photometry and nightly stacked photometry can help reveal the behavior of faint transients is shown in Figure 2. We note that public ZTF images acquired after Public Data Release 2 on 2019 December 11 were not yet released for forced photometry when this work was carried out.

Zoom In Zoom Out Reset image size

KNe are expected to fade rapidly, as AT2017gfo clearly confirmed. We automatically reject all those candidates that do not show any $\gt 2\sigma$ evolution for 6, 8, and 10 days from the first detection in g, r, and i band, respectively. Forced photometry and nightly stacking increased the measured duration of some candidates because they allowed us to recover additional faint data points. Thus we apply another filter on stacked light curves to reject those candidates with total time span $\gt 14$ days and with time lag between first and last detection larger than 10, 12, and 14 days in g, r, and i band, respectively.

Where possible, we use two different linear fits before and after the brightest point to measure the rise and fade rates. In each band, we required a minimum time baseline of 3 hr between first and last detection for the data points to be fit. If the brightest point of the light curve occurred within 3 hr of the first or the last detection, a single linear fit was performed. Otherwise, two different linear fits were used before and after the brightest point to measure both the rise and fade rates.

Thanks to the rich optical/infrared observations available, it is tempting to use AT2017gfo as a prototypical KN to set the fading (or rising) thresholds to photometrically select KN candidates based on their photometric evolution. However, AT2017gfo is unlikely to be representative for its class, as deep follow-up and joint afterglow and KN fits of short GRBs have demonstrated (Gompertz et al. 2018; Ascenzi et al. 2019; Rossi et al. 2020). Instead, we chose thresholds for the fading rate based on a grid of KN models constructed with the Monte Carlo radiative transfer code POSSIS (Bulla 2019). Specifically, we use a subset of the BNS grid presented in Dietrich et al. (2020) and the NSBH grid presented in Anand et al. (2020). Both grids use geometries that are axially symmetric and constructed by varying three parameters: the mass ejected on dynamical timescales ("dynamical ejecta," ${M}_{\mathrm{ej}}^{\mathrm{dyn}}$), the mass released from the merger remnant and debris disk ("post-merger" ejecta, ${M}_{\mathrm{ej}}^{\mathrm{pm}}$), and the observer viewing angle (${{\rm{\Theta }}}_{\mathrm{obs}}$, with ${{\rm{\Theta }}}_{\mathrm{obs}}=0^\circ$ corresponding to a face-on view of the system). KN light curves are predicted for 11 viewing angles from a polar (${{\rm{\Theta }}}_{\mathrm{obs}}=0^\circ$) to an equatorial (${{\rm{\Theta }}}_{\mathrm{obs}}=90^\circ$) orientation and for the following choices of the ejecta masses: ${M}_{\mathrm{ej},\mathrm{BNS}}^{\mathrm{dyn}}\,=[0.001,0.005,0.01,0.02]\,{M}_{\odot }$, ${M}_{\mathrm{ej},\mathrm{BNS}}^{\mathrm{pm}}\in [0.01,0.13]\,{M}_{\odot }$ with a mass step of $0.02\,{M}_{\odot }$, and ${M}_{\mathrm{ej},\mathrm{NSBH}}^{\mathrm{dyn}}$ and ${M}_{\mathrm{ej},\mathrm{NSBH}}^{\mathrm{pm}}\,\in [0.01,0.09]\,{M}_{\odot }$ with a mass step of $0.01\,{M}_{\odot }$. A half-opening angle of 30° is assumed for the lanthanide-rich dynamical component in both systems. This results in 308 and 891 KN light curves for the BNS and NSBH grids, respectively. We refer the reader to Dietrich et al. (2020) and Anand et al. (2020) for more details about the adopted BNS and NSBH grids.

In Figure 3, we show the distribution of fading rates calculated between light curve peak epoch and peak time +6 days. In this work, we considered the threshold that constrains 95% of models of the grid for each filter, i.e., 0.57, 0.39, and 0.27 mag day⁻¹ for g, r, and i band respectively. All those threshold choices fully encompass the fading rate of AT2017gfo (see blue vertical lines in Figure 3). Color evolution was not considered as a rejection criterion a priori, but it is discussed for the selected sources in Section 3.

Zoom In Zoom Out Reset image size

2.4. Candidate Vetting

In this section, we present the results of our searches, divided into an all-sky search and a galaxy-targeted search. A total of 24 sources survived the selection criteria and strict quality checks, but none of them was deemed a robust KN candidate.¹⁸ The main properties of the selected candidates are presented in Table 1, where candidates are grouped by (i) possible crossmatch with CLU galaxies, (ii) proximity to an extended source further than 15 and no underlying source, (iii) presence of an underlying source within 15, and (iv) "hostless" candidates. Appendix B presents the candidates' decay rates in each band and a detailed post mortem of individual candidates can be found in Appendix A.

Table 1.  Candidates that Passed Our Selection Criteria Described in Section 2.3

Name	R.A.	Decl.	${b}_{\mathrm{gal}}$	distnr	sgscore	CLU Match	Comment
Name	(deg)	(deg)	(deg)	(arcsec)		(Boolean)
ZTF18aazjzeda	245.80994	65.01023	39.8	⋯	⋯	1	Cataclysmic variable
ZTF19abzwbxyb	305.86048	6.66687	−17.0	<1.5	0.13	1	QSO
ZTF18achqdzhc	40.20850	−9.17101	−58.5	<1.5	0.12	1	Bright nuclear variability
ZTF20aaouvjnb	169.87812	15.41104	65.7	⋯	⋯	1	Multiple outbursts in ZTF
ZTF19abqtcoba	295.23436	2.22770	−10.0	⋯	⋯	1	Extragalactic nova or foreground cataclysmic variable?
ZTF18abyzkeq	342.39259	37.93156	−19.0	14.35	0.10	0	Cataclysmic variable
ZTF19aayhxlo	295.77816	48.32182	12.0	8.85	0.06	0	Post-detections in ZTF
ZTF19acecsfi	119.87837	44.27264	30.3	23.62	0.02	0	Likely stellar/cataclysmic variable
ZTF19aabgebm	29.88002	23.84546	−36.4	9.09	0.01	0	GRB 190106A afterglow
ZTF20aahcrjn	109.06754	35.75254	20.2	9.70	0.03	0	Multiple bursts in ZTF
ZTF18abuzpri	288.24987	61.55687	21.2	<1.5	0.32	0	Stellar
ZTF19abqiwjq	293.06491	6.72184	−5.9	<1.5	0.38	0	Pre-activity in PS1 images
ZTF19abqneae	295.60303	0.56979	−11.1	<1.5	0.00	0	Pre-activity in PS1 images
ZTF19acbtthv	16.02320	50.12632	−12.7	<1.5	0.49	0	Pre- and post-detections in PS1 and ZTF
ZTF19acszwgx	271.89500	10.59102	14.5	<1.5	0.03	0	Pre-activity in PS1 images
ZTF19abcputm	294.13549	−7.35038	−13.4	<1.5	0.14	0	Pre-activity in PS1 images
ZTF19abogdfr	293.39693	−17.64588	−17.1	<1.5	0.08	0	Pre-activity in PS1 images
ZTF19aapbbde	270.30443	13.28065	17.0	<1.5	0.37	0	Pre-activity in PS1 images
ZTF19abcpiag	264.64519	1.42066	16.8	⋯	⋯	0	Cataloged variable source
ZTF18acfmhrt	30.63920	55.41712	−6.1	⋯	⋯	0	Color evolution unlike KN
ZTF19abxwwmr	313.90428	22.98274	−14.2	⋯	⋯	0	Color evolution unlike KN
ZTF19aanhtzz	200.71300	57.46357	59.2	⋯	⋯	0	Cosmological afterglow?
ZTF19abudvoz	46.68777	46.15358	−10.6	⋯	⋯	0	Pre-activity in PS1 images
ZTF18acsjqjd	29.91694	54.17361	−7.4	⋯	⋯	0	Pre-activity in PS1 images

Notes. Along with the equatorial J2000 coordinates and Galactic latitude (${b}_{\mathrm{gal}}$), we indicate the distance to the closest source within $30^{\prime\prime}$ with star/galaxy classification score $\mathrm{sgscore}\leqslant 0.1$, suggesting the source to be a possible galaxy. The star/galaxy classification score of those candidates lying on top of faint sources ($\mathrm{distnr}\lt 1.5^{\prime\prime}$) is reported even where $\mathrm{sgscore}\gt 0.1$. The 24 candidates are grouped by their possible crossmatch with CLU galaxies, proximity to an extended source further than 15 and no underlying source, presence of an underlying source within 15, and finishing with "hostless" candidates. Finally, the comment column summarizes the post mortem description of each candidate from Section A.

^aWithin 100 kpc from a galaxy <10 Mpc away. ^bFound in both all-sky and galaxy-targeted searches. ^cFound only in galaxy-targeted search.

Download table as:  ASCIITypeset image

3.1. All-sky Search

The untargeted search in 514 fields, requiring $\geqslant 3$ detections, yielded 10,419 candidates, 23 of which passed our thorough vetting. Of these, four are spatially consistent with cataloged variable sources, 10 show variability in archival Pan-STARRS 1 (PS1; Chambers et al. 2016) data or in additional ZTF data, three are likely stellar in origin, and three were rejected mainly because of the combination of color and color evolution (Figure 4) being incompatible with our grid of KN models. All the candidates rejected because of their color evolution were found at relatively low Galactic latitude $| {b}_{\mathrm{gal}}| \leqslant 20$ deg, which supports the hypothesis of a stellar origin. This leaves us with three candidates.

Zoom In Zoom Out Reset image size

ZTF19abqtcob/AT2019aadh is located 55 kpc (2188) from a galaxy 8.6 Mpc away. The light curve (Figure 5) shows a rapid r-band decay of ${\alpha }_{r}=0.43$ mag day⁻¹, but no color information is available. Assuming an association with the nearby galaxy, the absolute magnitude at peak would be ${M}_{r}\sim -10$, consistent with classical nova luminosity. In addition, it is possible that ZTF19abqtcob is a foreground Galactic variable because of its low latitude (${b}_{\mathrm{gal}}=-10.0$) and the lack of any other possible host galaxy visible in its proximity. Therefore, we refrain from considering ZTF19abqtcob a viable KN candidate.

Zoom In Zoom Out Reset image size

The two remaining candidates, ZTF19aabgebm/AT2019aacx and ZTF19aanhtzz/AT2019aacu, are both located at high Galactic latitude ($| {b}_{\mathrm{gal}}| \gt 36$) and show red ($g-r\gt 0$ mag) color (Figure 5). ZTF19aabgebm is temporally and spatially coincident with the afterglow of the long GRB 190106A (Sonbas et al. 2019; Yurkov et al. 2019) at redshift z = 1.859 (Schady et al. 2019). Forced photometry allowed us to calculate a decay rate for ZTF19aabgebm of (1.24 ± 0.26)mag day⁻¹ in both g and r bands. An independent identification of ZTF19aabgebm in ZTF data is described in Ho et al. (2020a).

ZTF19aanhtzz could be considered the most puzzling transient found during our searches. The red color ($g-r\sim 0.3$ mag) disfavors the Galactic cataclysmic variable scenario, but there is no recorded gamma-ray signal temporally and spatially compatible with the optical transient to support the cosmological afterglow hypothesis. Deep g- and r-band images were acquired with the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) at W. M. Keck Observatory on 2020 June 22 UT, with 600 s of exposure time (Figure 6). LRIS data were processed using lpipe (Perley 2019), a fully automatic data reduction pipeline for imaging and spectroscopy. Photometry was performed using SExtractor (Bertin & Arnouts 2010) and it resulted in the non-detection of a source at the transient location to 3σ limiting magnitudes of $G\gt 26.8$ and $R\gt 27.1$. A source is present at an angular distance of 35 in PS1 and Legacy Survey DR8 images ($r\sim 22$, labeled as source "A" in Figure 6) and it can be modeled with a PSF shape according to the Legacy Survey DR8 catalog (Dey et al. 2019), which suggests the source to be stellar. No extended emission from source A is clearly detectable in the LRIS images.

Zoom In Zoom Out Reset image size

The galaxy WISE J132245.14+572830.4, labeled as source "B" in Figure 6, is placed at a distance of 267 Mpc by the Galaxy List for the Advanced Detector Era catalog (GLADE; Dálya et al. 2018), which would imply a projected distance of ∼83 kpc between the galaxy and ZTF19aanhtzz, if it was indeed the transient's host. Such a distance is within the range of projected offsets for neutron star mergers, considering short GRBs as a proxy (Berger 2014). On 2020 July 2 UT, we acquired a spectrum of WISE J132245.14+572830.4 using the Double Spectrograph (DBSP) mounted on the Palomar 200 inch telescope with 600 s of exposure time. The data were processed with a custom PyRAF reduction pipeline for DBSP (Bellm & Sesar 2016). The spectrum shows prominent Balmer series and other common galaxy features that place the galaxy at redshift z = 0.103, i.e., at a luminosity distance of 490.57 Mpc and a projected distance of ∼135 kpc from the transient. The projected distance is too large for WISE J132245.14+572830.4 to be a reliable host of ZTF19aanhtzz (although short GRBs were found with a projected distance from the host larger than 100 kpc; Berger 2014).

If neither source A nor source B are the transient's host, we consider the hypothesis of an underlying host galaxy being too faint to be detectable in our Keck+LRIS images. If such a faint host was located at a distance of <250 Mpc, its luminosity would be $M\sim -10$ or fainter, typical of dwarf galaxies satellite of more massive galaxies. The apparent lack of such a massive galaxy in the field of ZTF19aanhtzz makes this scenario improbable. It is possible that a more massive host galaxy at high redshift exists that was not detectable in our Keck+LRIS images. If ZTF19aanhtzz was a Galactic source, located in the outskirts of the Milky Way, the LRIS limit suggests the absolute magnitude of the quiescent counterpart to be $M\gt 7$ mag. Based on a typical HR diagram, this only allows a main-sequence star that is later than K0 in spectral type, or a white dwarf, to be the progenitor star, but an M-dwarf progenitor is likely excluded by the non-detection in z band by the Legacy Survey. We consider now the scenario in which ZTF19aanhtzz was a cosmological relativistic afterglow. The lack of a detected GRB and the observation of the rise phase of the transient make it unlikely for ZTF19aanhtzz to be the afterglow of a relativistic explosion viewed on-axis. Although a detail analysis using off-axis GRB models (e.g., van Eerten et al. 2010; Beniamini et al. 2020) is outside the scope of this paper, we suggest that a slightly off-axis relativistic explosion might explain the behavior of this transient. In summary, with the data in hand, a conclusive answer regarding the nature of ZTF19aanhtzz/AT2019aacu is yet to be found, but a cosmological afterglow origin scenario is favored.

3.2. Galaxy-targeted Search

In this work, we used rigorous criteria to identify uncommon transients that could possibly be KNe. The rich ZTF data set provided us with high cadence photometry and, usually, color information that we could use to find rapidly fading candidates possibly associated with nearby galaxies and separate them from other classes of transients or variables. Our searches did not yield any source that we claim to be a viable KN candidate, although the nature of some fast transients remains unknown, including the possibly extragalactic fast transient ZTF19aanhtzz. It is possible that KNe were present in the data set, but they were missed in this analysis because, for example, they were too close to the detection limit, due to the pipeline detection efficiency, or because they evolved in different ways than we expected based on available models and on the knowledge of GW170817.

This search for KNe in the ZTF survey, independent of GW or short GRB triggers, complements the ZTF searches for optical counterparts to 13 neutron start mergers found in GWs during O3, described in Kasliwal et al. (2020). In that work, we used the query and photometry methods described in Section 2 to search for KNe by requiring $\geqslant 1$ detection in the alert stream over the 95% of the localization skymaps in the 3 days following the mergers. The analysis presented in Kasliwal et al. (2020) led to no viable KN counterpart to O3 neutron star mergers.

This work further demonstrated the importance of repeated, deep observations of wide-field surveys to better understand the nature of newly identified astronomical transients. A key problem in KN searches in "archival" survey data is the ability to recognize KNe without spectroscopic follow-up. The problem of photometric identification of KNe among other transients was already addressed, for example, by Cowperthwaite & Berger (2015), Doctor et al. (2017), Cowperthwaite et al. (2018), and Bianco et al. (2019). In general, these works find that a combination of color information (at a given time) and rapid evolution is effective at separating KNe from other extragalactic sources such as supernovae. In our case, the selection criteria on the fading rate and transient duration had already helped us remove all supernova candidates from our sample, leaving AGN and Galactic sources as main contaminants, as expected from previous fast transient searches (e.g., van Roestel et al. 2019; Andreoni et al. 2020a).

We use the grid of models described in Section 2.3 to visualize the expected g − r color evolution over time for BNS and NSBH mergers. Magnitudes and colors (specifically g − r in this case) are expressed in the AB system. Figure 4 shows that all those candidates for which inter-night color evolution is available have color evolution different than predicted by KN models. Interestingly, the candidate that approaches KN models the most is ZTF19aabgebm, the red and rapidly evolving afterglow of GRB 190106A. We also considered a set of 30 KN models from the grid computed by Kasen et al. (2017) with parameters similar to the blue component of AT2017gfo in the ranges ${M}_{\mathrm{ej}}\in [0.01,0.04]\,{M}_{\odot }$, ${v}_{k}\in [0.1\,-0.3]$, and ${X}_{\mathrm{lan}}\in [1{\rm{e}}-5,1{\rm{e}}-4]$ for ejecta mass, kinetic ejecta velocity, and lanthanide fraction, respectively. In the first 3 days from the merger, these models are characterized by blue colors and thus provide better match to our KN candidates compared to the POSSIS grid shown in Figure 4. However, the resulting light curves get redder even more steeply than the POSSIS grid, placing the candidates even further from the region of the phase space expected for KNe at later epochs.

In general, systematic searches in the broadest possible parameter space covering the models presented in Section 2.3 yield reasonable fits to our selected candidates (Coughlin et al. 2018; Stachie et al. 2019); the lack of either explosion time or distance constraints make it possible to fit them within error bars of the models. This emphasizes the need for rapid photometric follow-up of sources of this type to make the most straightforward classifications.

The searches conducted during the first 23 months of the ZTF survey, based on alerts, have not yielded any likely KN candidate. This result can be used to constrain KN rates to compare against neutron star merger rates obtained from GW and short GRB observations.

Rates were estimated using simsurvey¹⁹ (Feindt et al. 2019), a simulation software package that accounts for custom transient light curves and observational parameters of the survey. simsurvey uses Monte Carlo methods to simulate synthetic light curves. We used ZTF data to calculate the number of KNe detected for each rate used in the injection of synthetic transients. Milky Way extinction was applied based on the Schlegel et al. (1998) reddening maps. We considered three types of KN light curves: a 3 days constant luminosity (or "top-hat") model, a linear decay model, and models based on radiative transport simulations. Synthetic sources are injected in simsurvey assuming a broad range of rates. Our limits on the rate are chosen at the point where three of the injected KNe are recovered, as dictated by Poissonian small number statistics for a 95% confidence level (Gehrels 1986).

• 1.  
  Top-hat—Constant 3 days light curve model (Appendix C). The top-hat model at ${M}_{g,r,i}=-16$ (about the maximum brightness measured for AT2017gfo) yielded a rate of $R\lt 398$ Gpc⁻³ yr⁻¹ for at least one 5σ detection, $R\lt 930$ Gpc⁻³ yr⁻¹ for at least two 5σ detections separated by at least 3 hr, and $R\lt 1121$ Gpc⁻³ yr⁻¹ by requiring at least three 5σ detections. Appendix C presents a table with the results for a range of absolute magnitudes $M\in [-12,-17]$.
• 2.  
  Linear decay—We considered a grid of linear decay models by varying the starting absolute magnitude $M\in [-14.5,-17.5]$ mag and the decay rate $\alpha \,\in [0.3,1.5]$ mag day⁻¹ to encompass the most probable decay rates that we considered to chose our selection criteria (Figure 3, Appendix D). These models do not account for color evolution and are filter-agnostic. The table in Appendix D presents the results for the linear decay grid and Figure 7 represents contours indicating rates of R = 500, 1000, and 2000 Gpc⁻³ yr⁻¹.
• 3.  
  Radiative transport simulations—We use the best-fit model for AT2017gfo obtained in Dietrich et al. (2020) and assume two KN populations: the first where all KNe have light curves like AT2017gfo, the second where KNe are instrincally like AT2017gfo, but with viewing angle i uniformly distributed in cos(i). The results are shown in Appendix E. In addition, we explore the parameter space of ejecta masses for the models tailored to BNS and NSBH from Dietrich et al. (2020) and Anand et al. (2020), also described in Section 2.3. The resulting rates for the grid of models obtained with a fixed half-opening angle ϕ = 30 deg are shown in Appendices F and G and illustrated in Figure 8. We note that the NSBH grid provides less constraining results. The light curves are fainter at peak, and broader compare to the BNS grid. Therefore, the detected light curves would be observable for longer times, but in a smaller distance range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The rates obtained by requiring three 5σ detections represent our all-sky search conservatively. In fact, we queried the ZTF database selecting sources with at least three detections in the alerts. However, only one alert (chronologically the last one) needs to have $\geqslant 5\sigma$ significance. For example, it is possible that only one alert is issued for a transient when its flux passes the 5σ threshold, but the alert packet includes information about low-significance detections (between 3σ and 5σ) that did not trigger any alert generation. These low-significance detections were accounted for during our query. In addition, our experiment targeting CLU galaxies required only two detections at least 3 hr apart, the first of which could be a low-significant detection below the 5σ mark.

On the other hand, simsurvey does not account for the detection efficiency curve of the ZTF data processing pipeline, which is currently unavailable. Frohmaier et al. (2017) calculated the detection efficiencies for the Palomar Transient Factory (PTF) pipeline (that relies on the image-subtraction algorithm developed by Alard & Lupton 1998), obtaining a completeness of $\sim 95 \%$ where the host galaxy surface brightness is ≲10% of the transient brightness. From the PTF experience, we can also expect the detection completeness to drop to 60%–70% in the presence of bright hosts with $r\lesssim 18$ mag (Frohmaier et al. 2017). However, a dedicated study of the recovery efficiency of the ZTF pipeline (that relies on the ZOGY image-subtraction algorithm; Zackay et al. 2016) is necessary to provide exact completeness figures. Such a study lies outside the scope of this paper.

We estimated that the random misalignment between science and reference images can cause a loss of ∼1% of the chip area. Also, Fremling et al. (2020) calculated that <1% of astronomical transients will be lost due to proximity to bright sources outside the Galactic plane.

Catalog incompleteness should also be taken into account. A galaxy number completeness lower limit of 57% integrated between 0 < z < 0.05 was obtained for the NASA/IPAC Extragalactic Database (NED) by Fremling et al. (2020) using Type Ia supernovae as proxy. The CLU catalog encompasses NED and includes >18,000 galaxies within 200 Mpc found during the PTF Hα survey (Cook et al. 2017), so we expect our integrated completeness to be larger than 57%.

With these caveats in mind, we will now discuss the main results obtained in this analysis, considering the case of two 5σ detections as the best representation of our results.

Using the best-fit model to AT2017gfo in radiative transfer simulations, for two 5σ detections, we obtained a rate of R < 1775 Gpc⁻³ yr⁻¹, while considering a population of KNe with the same intrinsic properties of AT2017gfo and with a uniform viewing angle distribution, we obtained a rate of R < 4029 Gpc⁻³ yr⁻¹. As discussed in Section 3, models obtained with POSSIS radiative transport simulations tend to err toward redder and fainter KNe. Assuming a blue KN with a linear decay rate of α = 0.7 mag day⁻¹ as for AT2017gfo between g and r band with a peak magnitude of M = −16.5, we can constrain the KN rate to be $R\lt 937$ Gpc⁻³ yr⁻¹ (1363 Gpc⁻³ yr⁻¹) by requiring two (three) detections in ZTF data (5σ).

Before the discovery of GW170817, Doctor et al. (2017) constrained the bright KN rate to be $R\lt 2.4\times {10}^{4}$ Gpc⁻³ yr⁻¹ (90% confidence) with the Dark Energy Survey. Then, other surveys constrained the KN rate using AT2017gfo as reference. Smartt et al. (2017) set the KN rate to be $R\lt 3.0\times {10}^{4}$Gpc⁻³ yr⁻¹ (95% confidence) with ATLAS and Yang et al. (2017) placed an upper limit of $R\lt 9.9\times {10}^{4}$ Gpc⁻³ yr⁻¹ (90% confidence) with 13 months of nightly monitoring of nearby galaxies with DLT40. The limits calculated in this work are more constraining than those listed above. A more stringent limit of $R\lt 800$ Gpc⁻³ yr⁻¹ was placed by Kasliwal et al. (2017) with PTF data, also focusing on galaxies present in the CLU catalog (although a correction for the galaxy catalog completeness was not applied). Our results are consistent with ZTF survey simulation predictions described in Sagués Carracedo et al. (2020). The non-detection of a KN in almost 2 yr of survey suggests that the number of KNe detectable with ZTF estimated by Scolnic et al. (2018) (using a fixed rate of R = 1000 Gpc⁻³ yr⁻¹ and average values for the cadence, zero-point, and sky noise) may have been too optimistic.

When comparing KN rates with BNS merger or short GRB rates, a KN luminosity function should be considered, but its uncertainties are still large (Kasliwal et al. 2020). Our limits using a GW170817-like model are consistent with isotropic equivalent short GRB rates including, for example, $R={270}_{-180}^{+1580}$ Gpc⁻³ yr⁻¹ (Fong et al. 2015), $R={352}_{-281}^{+810}$ Gpc⁻³ yr⁻¹ (della Valle et al. 2018), or $R={160}_{-100}^{+200}$ Gpc⁻³ yr⁻¹ (Dichiara et al. 2020). ZTF results are also consistent with rates estimated with population synthesis models such as $R={600}_{-300}^{+600}$ Gpc⁻³ yr⁻¹ (Chruslinska et al. 2018), a solution among those that returned higher rates that are consistent with the BNS merger rate from GWs (Abbott et al. 2017b, 2020) (see also Tutukov & Yungelson 1993).

Assuming that all BNS mergers generate a bright optical transient, the ZTF limits can constrain the upper side of rates obtained via GW observations of 250–2810 Gpc⁻³ yr⁻¹ (Abbott et al. 2020), under the assumption that most KNe are intrinsically similar to AT2017gfo. The ZTF KN rate can also constrain the rates of $R={1109}_{-657}^{+1432}$ Gpc⁻³ yr⁻¹ (Jin et al. 2018), calculated from short GRB observations, and $R={830}_{-661}^{+2091}$ Gpc⁻³ yr⁻¹, derived from BNS system observations in the Milky Way (${R}_{\mathrm{BNS}}={83.0}_{-66.1}^{+209.1}$ Myr⁻¹ at 95% confidence level; Kalogera et al. 2004) assuming a galaxy number density of $\sim {10}^{-2}$ Mpc⁻³ as in della Valle et al. (2018). Our results are consistent with BNS merger rates obtained by Kim et al. (2015), who calculated ${R}_{\mathrm{BNS}}={21}_{-14}^{+28}$ Myr⁻¹ at 95% confidence level. Similarly, our result is compatible with the updated Milky Way BNS merger rate obtained by Pol et al. (2019, 2020) of ${R}_{\mathrm{BNS}}={37}_{-11}^{+24}$ Myr⁻¹ at 90% confidence level. The rates from Kim et al. (2015) and Pol et al. (2020) translate into $R={210}_{-140}^{+280}$ Gpc⁻³ yr⁻¹ and $R={370}_{-110}^{+240}$ Gpc⁻³ yr⁻¹, respectively, using again a galaxy number density of 10⁻² Mpc⁻³.

Figure 9 offers a visual representation²⁰ of the rates discussed in this section.

Zoom In Zoom Out Reset image size

If changes in the viewing angle make the KN rate become R < 4029 Gpc⁻³ yr⁻¹ as we predicted using radiative transfer models (Bulla 2019), our ZTF limits would become unconstraining for BNS mergers. The current rate from GWs does not account for NSBH mergers, which also are expected to be accompanied by bright KNe under favorable conditions of mass ratios and neutron star radii, thus it could be interesting to compare the ZTF constrains with new BNS + NSBH merger rates from GW observations in the near future.

We explored 23 months of ZTF survey data searching for rapidly fading transients that could be classified as KNe based on their photometric evolution. In addition to ZTF alert database queries, we employed forced PSF photometry and nightly flux stacking to improve our ability of recognizing KNe among $\gt {10}^{4}$ candidates. All most promising transients turned out to be either Galactic foreground or cosmological background sources, leaving no viable KN in the sample. The observations allowed us to place constraints on the KN rate for a variety of radiative transport-based models, specifically $R\lt 1775$ Gpc⁻³ yr⁻¹ for KNe similar to AT2017gfo and R < 4029 Gpc⁻³ yr⁻¹ for KNe with the same intrinsic properties of AT2017gfo assuming a uniform viewing angle distribution.

Several ongoing surveys observe large volumes of sky at increasingly high cadence including, for example, Pan-STARRS (Chambers et al. 2016), the Asteroid Terrestrial-impact Last Alert System (ATLAS²¹ ), and All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Holoien et al. 2017), ZTF itself, and multi-facility programs such as "Deeper, Wider, Faster" (Andreoni et al. 2020a; J. Cooke et al. 2020, in preparation). Dedicated searches in these large data sets may be necessary to find KNe (see for example Smartt et al. 2019; McBrien et al. 2020, regarding ongoing KN searches with PS1).

Based on this study, it is clear that both increased depth and rapid follow-up of interesting candidates, as executed for example in the ZTF follow-up of GW events (Coughlin et al. 2019a; Kasliwal et al. 2020, and Anand et al. 2020), should be prioritized to maximize the probability of a KN discovery. In about 2 yr, ZTF has allowed us to place tighter constraints on the KN rate than most other surveys, but there is room for improvement in the near future. Two more years of ZTF survey will naturally halve the rate limits that we calculated. To maximize the volume explored, depth should be prioritized over area (e.g., Ghosh et al. 2017). Nightly stacking of ZTF images (D. A. Goldstein et al. 2020, in preparation) or longer exposures may be successful avenues for KN discovery. The identification of KN candidates in real time is indisputably valuable because it can trigger rapid multi-wavelength follow-up observations, including optical/infrared spectroscopy. On the other hand, the limited availability of 8 m class telescopes, necessary to classify sources fainter than ∼21 mag, represents a challenge for discovery programs relying on deeper observations.

This work has shown that photometric monitoring combined with archival information from multiple surveys can also be very effective at separating promising KN candidates from other types of rapidly evolving transients. The ability to recognize KNe solely based on their light curves will be key during the upcoming Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) at the Vera C. Rubin Observatory. Photometric KN candidate identification will be possible, in particular, if a "rolling" cadence is chosen that entails nightly multi-band imaging (Andreoni et al. 2019a). Although ongoing surveys have some potential to detect KNe independently of GW or GRB triggers, we expect high-cadenced LSST observations to yield many serendipitous KN discoveries.

We thank the anonymous referee who helped us improve the quality of the manuscript. We thank P. Nugent for reviewing the manuscript and providing us with useful feedback. We thank D. Perley for useful discussion. Based on observations obtained with the Samuel Oschin 48 inch Telescope and the 60 inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project, a scientific collaboration among the California Institute of Technology, the Oskar Klein Centre, the Weizmann Institute of Science, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron, the University of Wisconsin-Milwaukee, and the TANGO Program of the University System of Taiwan. Further support is provided by the U.S. National Science Foundation under grant No. AST-1440341.

This work was supported by the GROWTH (Global Relay of Observatories Watching Transients Happen) project funded by the National Science Foundation under PIRE grant No. 1545949. GROWTH is a collaborative project among California Institute of Technology (USA), University of Maryland College Park (USA), University of Wisconsin-Milwaukee (USA), Texas Tech University (USA), San Diego State University (USA), University of Washington (USA), Los Alamos National Laboratory (USA), Tokyo Institute of Technology (Japan), National Central University (Taiwan), Indian Institute of Astrophysics (India), Indian Institute of Technology Bombay (India), Weizmann Institute of Science (Israel), The Oskar Klein Centre at Stockholm University (Sweden), Humboldt University (Germany), Liverpool John Moores University (UK) and University of Sydney (Australia).

E.C.K., A.S.C., and A.G. acknowledge support from the G.R.E.A.T research environment funded by Vetenskapsrådet, the Swedish Research Council, project number 2016-06012. E.C.K. also acknowledges support from The Wenner-Gren Foundations. M.W.C. acknowledges support from the National Science Foundation with grant No. PHY-2010970. D.L.K. is supported by NSF grant AST-1816492. M.M.K. acknowledges generous support from the David and Lucille Packard Foundation.

We thank A. Miller for the star/galaxy classification code development for ZTF. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (doi:10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Software: Astropy (The Astropy Collaboration 2018), Pandas (pandas development team 2020), Matplotlib (Caswell et al. 2020), Kowalski (Duev et al. 2019), lpipe (Perley 2019), PostgreSQL²² , PyRAF DBSP data reduction pipeline (Bellm & Sesar 2016), ztfquery (Rigault 2018), SExtractor (Bertin & Arnouts 2010), ForcePhotZTF (Yao et al. 2019).

We discuss here individual candidates. The main features that we consider are: (1) variable catalog crossmatches; (2) possible CLU galaxy crossmatch, or presence of a bright nearby galaxy; (3) past variability in survey data other than ZTF; (4) color evolution.

• 1.  
  ZTF18aazjzed—Reported to Transient Name Server (TNS) as AT2018ltl, it was cataloged as a cataclysmic variable. There are no PS1 pre-detections, but a faint source is visible in the ZTF reference image and there are recent ZTF detections (after the end date of the experiment). In addition, the color evolution is incompatible with KN models (Figure 4).
• 2.  
  ZTF19abzwbxy—The transient is located 86.5 kpc from a CLU galaxy 67.2 Mpc away. However, it is coincindent with a faint source present in the PS1 catalog, making an association with the CLU galaxy unlikely. In addition, the source is cataloged as variable in the AAVSO International Variable Star Index, possibly being a quasi-stellar object (QSO; Brescia et al. 2015).
• 3.  
  ZTF18achqdzh—Nuclear source located in a galaxy with spectroscopic redshift z = 0.0461, so at peak the transient is at ${M}_{r}\lt -17$, probably too bright to be a KN (Kasliwal et al. 2020). The host was classified as a low-ionization nuclear emission-line region (LINER) and can harbor an AGN (Cavuoti et al. 2014). Sporadic faint detections in ZTF forced photometry suggest the emission to be due to nuclear variability.
• 4.  
  ZTF20aaouvjn—Reported to TNS as AT2020lru. Located 76.2 kpc (388'') from a galaxy 40.5 Mpc away. About 25 from a star, most likely to be the source if it has large proper motion. No other sources are present nearby and there are no PS1 pre-detections. The color evolution is incompatible with KN models (Figure 4).
• 5.  
  ZTF19abqtcob—Reported to TNS as AT2019aadh, it is located 55 kpc (1313'') from a CLU galaxy 8.6 Mpc away. If the host is indeed the nearby CLU galaxy, this would imply an absolute magnitude of $M\sim -10$ at peak, consistent with a classical nova. The low Galactic latitude of ${b}_{\mathrm{gal}}=-10.0$, and the lack of any other possible host nearby make it likely for ZTF19abqtcob to be a foreground Galactic cataclysmic variable.
• 6.  
  ZTF18abyzkeq—Hostless source with no PS1 pre-detections, ZTF18abyzkeq is classified as a variable star in the AAVSO International Variable Star Index. As expected, the color evolution is incompatible with KN models (Figure 4).
• 7.  
  ZTF19aayhxlo—Reported to TNS as AT2019aada, ZTF19aayhxlo appears to be a hostless source with no PS1 pre-detections. The color evolution is incompatible with KN models (Figure 4). Recent ZTF detections (after the end date of the experiment) suggests the source to be a variable object.
• 8.  
  ZTF19acecsfi—Faint PS1 and SDSS pre-detection within $1^{\prime\prime}$. Deep Legacy Survey images suggest the presence of an underlying blue source deemed as stellar in Legacy Survey DR8 catalog (Dey et al. 2019). Despite the high Galactic latitude ${b}_{\mathrm{gal}}=30.3$, we consider the source to be likely a cataclysmic variable.
• 9.  
  ZTF19aabgebm—Reported to TNS as AT2019aacx, ZTF19aabgebm is the (already known) optical afterglow of long GRB 190106A (Sonbas et al. 2019; Yurkov et al. 2019) at redshift z = 1.859 (Schady et al. 2019). We note that its color and color evolution are barely compatible with KN models for NSBH mergers (Figure 4).
• 10.  
  ZTF20aahcrjn— Reported to TNS as AT2020atp. A faint source is present at the candidate location in the ZTF reference image. There are no PS1 pre-detections, but forced photometry revealed another possible outburst in 2018 April, suggesting this to be a Galactic variable. Moreover, the color evolution is incompatible with KN models (Figure 4).
• 11.  
  ZTF18abuzpri—Reported to TNS as AT2018ltn. A prominent source is present in the ZTF images, likely consituted of two stars (Legacy Survey DR8). Although there is no activity recorded in the PS1 catalog, the likely stellarity and the significantly blue color $g-r=-0.73\pm 0.18$ makes it unlikely for ZTF18abuzpri to be a KN.
• 12.  
  ZTF19abqiwjq—Reported to TNS as AT2019aade, it is located in a crowded field at low Galactic latitude ${b}_{\mathrm{Gal}}=-5.9$ deg. A faint source is present in the reference image at the location of the candidate. PS1 detected this source multiple times in different filters and the flux in these detections has increased significantly (more than eight times), indicating past activity. We therefore rule out ZTF19abqiwjq as a KN candidate.
• 13.  
  ZTF19abqneae—Reported to TNS as AT2019aadd, it is located in a crowded stellar field at ${b}_{\mathrm{Gal}}=-11.1$ deg. The candidate has multiple PS1 detections on different filters and the flux in these detections has doubled, showing activity in the past.
• 14.  
  ZTF19acbtthv—Reported to TNS as AT2019aacw. Hostless source, however PS1 pre-detections are present within $2^{\prime\prime}$ and recent additional ZTF detections after the nominal end of the experiment, indicating repeated activity. Moreover, the color evolution is incompatible with KN models (Figure 4).
• 15.  
  ZTF19acszwgx—Previously reported to TNS by ATLAS team as AT2016hnn. An underlying faint source is visible in the ZTF reference image, with multiple PS1 detections within $1^{\prime\prime}$. The g-band PS1 flux varies 30% thus it was active in the past and a possible re-brightening could be detected in recent ZTF data.
• 16.  
  ZTF19abcputm—Reported to TNS as AT2019aacy. A faint source is present in the ZTF reference image, and it has multiple PS1 detection within $1^{\prime\prime}$. The r-band PS1 flux varies 300% thus it was active in the past.
• 17.  
  ZTF19abogdfr—Reported to TNS as AT2019aacz, the source shows multiple PS1 detections within $1^{\prime\prime}$. The g-band PS1 flux varies by a factor of 80, thus it was active in the past.
• 18.  
  ZTF19aapbbde—The source shows an underlying blue point-like source. ZTF19aapbbde has numerous PS1 detections within $1^{\prime\prime}$, with flux measurements suggesting $\sim 50 \%$ increase in its brightness. This is grater than what would be expected within the PSF flux 3σ error, suggesting previous activity. The source shows dramatic reddening at a 3–4 days phase, which is incompatible with the expectations from KN models (Figure 4).
• 19.  
  ZTF19abcpiag—Previously reported by the Gaia team to TNS as AT2017dfl, this source is a cataloged variable star listed in the AAVSO International Variable Star Index. The transient presents an extremely rapid r-band light curve, but no color information is available.
• 20.  
  ZTF18acfmhrt—Reported to TNS as AT2018ltm, the source appears to be hostless with no pre-detections in PS1. An underlying faint blue source is visible in PS1 images. The color evolution is incompatible with KN models (Figure 4).
• 21.  
  ZTF19abxwwmr—Hostless and with no PS1 pre-detections. However, the color evolution is incompatible with KN models (Figure 4).
• 22.  
  ZTF19aanhtzz—Reported to TNS as AT2019aacu. The nearest source is 35 away and it is probably stellar (SDSS, Legacy Survey DR8). This source presents a number of interesting features, including red color, rapid evolution, and high Galactic latitude of ${b}_{\mathrm{gal}}=59.2$ deg. The lack of an apparent host galaxy suggests the source to be either a luminous cosmological transient (the afterglow of a GRB that went undetected by gamma-ray telescopes?) or some type of Galactic variable, rather than a kilonova. However, the red color disfavors the Galactic scenario and a conclusive answer regarding the nature of ZTF19aanhtzz/AT2019aacu is yet to be found.
• 23.  
  ZTF19abudvoz—Previously reported to TNS by the ATLAS team as AT2016ayj in 2016. Apparently orphan, it has two PS1 pre-detections within $1^{\prime\prime}$. In addition, the light curve is also unlike what is expected for KNe, with no color evolution during the rise time and with a g-band decay of ∼1 mag day⁻¹.
• 24.  
  ZTF18acsjqjd—The source shows pre-activity in PS1 images and therefore we exclude that it can be a KN.

In Table B1 we present the fading rates (mag/day) for the candidates that passed our selection criteria and quality checks. Specifically, we show for each band the index obtained by performing a linear fit from the brightest points light curves built with alerts, forced photometry, and nightly stacked forced photometry.

Table C1 presents kilonova rate upper limits for 3-day constant luminosities light curves corresponding to one, two, and three detections (5σ) in ZTF at 95% confidence.

Table D1 presents results for a grid of linear decay models used to inject synthetic light curves in simsurvey. The columns correspond to the starting absolute magnitude, the decay rate, and the rates corresponding to one, two, and three detections (5σ) in ZTF at 95% confidence.

Table E1 presents kilonova rate upper limits for the best-fit model to AT2017gfo (Dietrich et al., 2020) corresponding to one, two, and three detections (5σ) in ZTF at 95% confidence.

Table F1 presents kilonova rate upper limits in the parameter space of ejecta masses using a model grid tailored for BNS mergers with a fixed opening angle of ϕ = 30 deg (see Section 2.3 and Dietrich et al., 2020; Anand et al., in press) corresponding to one, two, and three detections (5σ) in ZTF at 95% confidence.

Table G1 presents kilonova rate upper limits in the parameter space of ejecta masses using a model grid tailored for NSBH mergers with a fixed opening angle of ϕ = 30 deg (see 2.3 and Dietrich et al., 2020; Anand et al., in press) corresponding to one, two, and three detections (5σ) in ZTF at 95% confidence.