An Empirical Limit on the Kilonova Rate from the DLT40 One Day Cadence Supernova Survey

Binary neutron star (BNS) systems (Hulse & Taylor 1975) have been studied with great interest by the astronomical community because of their connection with many open problems of astrophysics, from short gamma-ray bursts (GRBs) to r-process element production and from the physics of very high density matter to gravitational waves. The number of known BNSs today is limited to a dozen systems (Lattimer 2012) and the rate of BNS coalescences is known with orders of magnitude of uncertainty (Abadie et al. 2010; Coward et al. 2012; Berger et al. 2013; Petrillo et al. 2013; Drout et al. 2014;Siellez et al. 2014; de Mink & Belczynski 2015; Dominik et al. 2015; Fong et al. 2015; Jin et al. 2015; Kim et al. 2015; Vangioni et al. 2016; Abbott et al. 2017a; Scolnic et al. 2017). The rate of BNS coalescences can be constrained from the population of galactic pulsars (Kalogera et al. 2004), from modeling the evolution of the binary system (de Mink & Belczynski 2015; Dominik et al. 2015), from the cosmic abundance of r-process elements (Vangioni et al. 2016), or from measuring the rate of short GRBs, most likely produced at the moment of the coalescence (Berger et al. 2013).

An alternative method to constrain the rate of BNS mergers is to constrain the rate of kilonovae detected in supernova (SN) search surveys. Kilonovae are thought to be the ubiquitous, isotropically emitting counterparts to neutron star mergers. They are expected to eject at very high velocity a small mass (0.01–0.05 M_Sun) mainly made of high opacity r-process heavy elements; hence, they are predicted to have a "red" spectrum, to be faint at maximum light (M_V ∼ −16 mag), and to decline quickly over the course of 1–2 weeks (e.g., Metzger et al. 2010; Kasen et al. 2013; Piran et al. 2013). This is in contrast to the most common SNe, which evolve on 0–100 day timescales (e.g., see Figure 1.1 of Kasliwal et al. 2011).

One clear hindrance to calculating the kilonova rate has been the lack of observed kilonova events (with possible exceptions; see Jin et al. 2016 for a compilation) in ongoing transients searches. Some progress has been made by using theoretical kilonova light curves to calculate the upper limit of their rate in programs like the Dark Energy Survey (Doctor et al. 2017). Additionally, there have been several recent attempts to estimate the rate of fast optical transients that evolve on sub-day timescales (τ ∼ 0.5 hr to 1 day; Berger et al. 2013) all the way up to ∼10 day timescales (Drout et al. 2014).

With the discovery of AT 2017gfo (Abbott et al. 2017a, 2017b), we can directly constrain the rate of kilonovae by using its light curve as a template while allowing for possible diversity in kilonova light curves and a range of extinction values. Here, we present the rate estimate for kilonovae using our observed light curve of AT 2017gfo and the data from the ongoing Distance Less Than 40 Mpc (DLT40) supernova search (Tartaglia et al. 2017). DLT40 is an SN search that points to galaxies within D ≲ 40 Mpc with a one day cadence.⁶ Given the magnitude limit of the program (r ∼ 19 mag) DLT40 is well suited to detect a nearby kilonova event. An advantage of DLT40 is that we can directly use the light curve of AT 2017gfo obtained with the same instrumental setup (Valenti et al. 2017) to get a direct limit for similar transients in the DLT40 program.

This Letter is organized as follows. In Section 2, we describe the DLT40 search and the survey operation during the O2 run. In Section 3, we will highlight the steps necessary to measure the rates, while in Section 4 we will compare our results with previous rate estimates and discuss the future prospects on how the improve the rates with or independently from the next LIGO/Virgo observing run.

DLT40 is a one day cadence⁷ search for SNe, targeting galaxies in the nearby universe (D ≲ 40 Mpc) and designed for the goal of discovering ∼10 SNe per year within one day of explosion (Tartaglia et al. 2017). Running continuously since 2016 September, DLT40 observes ∼400–600 galaxies every night using a 0.41 m PROMPT telescope (Reichart et al. 2005) at the Cerro Tololo Inter-American Observatory (CTIO).

With a field of view (FoV) of 10 × 10 arcmin², DLT40 is suitable for mapping nearby galaxies down to a limiting magnitude of r ∼ 19 mag in 45 s unfiltered exposures. The DLT40 catalog includes galaxies from the Gravitational Wave Galaxy Catalog (GWGC; White et al. 2011) at a declination <20°, with an absolute magnitude M_B < −18 mag, galactic extinction A_V < 0.5 mag, and recessional velocity v_r <3000 km s⁻¹ (corresponding to D ≲ 40 Mpc). An additional constraint is that the FoV should not include stars brighter than 8 mag. These selections led to a sample of 2220 bright galaxies. In Figure 1, we show the sky distribution of galaxies in the DLT40 catalog. By comparing the integrated luminosity of DLT40's galaxy sample with the total integrated luminosity of the GWGC catalog, we estimated that the DLT40 catalog contains ∼60% of the total luminosity (and therefore roughly mass) of the GWGC catalog. Note that the GWGC is supposed to be complete out to 40 Mpc (White et al. 2011). The DLT40 galaxy catalog is not biased with respect to morphological type of galaxy (S. Yang et al. 2017, in preparation), but by construction it is biased against low-luminosity galaxies (for example, dwarf galaxies).

Zoom In Zoom Out Reset image size

In order to find SNe within a day from explosion, the DLT40 data are processed in nearly real time. Pre-reduced images are transferred within ∼1 minute to our dedicated server where new candidates are detected by image subtraction with respect to a template image (using Hotpants; Becker 2015) so that new candidates are available for scanning by human operator within a few minutes after acquisition. Confirmation images of new SN candidates are taken within a few hours of their DLT40 discovery, often with the Las Cumbres Observatory telescope network (Brown et al. 2013).

Since the beginning of the search (2016 September), we have discovered 26 SNe, 12 of which were first announced by DLT40. Seven were discovered within 48 h of explosion (see Table 1 for a list of transients discovered by DLT40). The late discovery of the remaining transients by DLT40 was due to poor weather conditions.

Table 2.  DLT40 Rate Estimation Results

Type	Extinction	Control Time	Lums Ratea	Vol Rateb	Milky Way Ratec
	(mag)	(days)	(SNuB)	(10−4 Mpc−3 yr−1)	((100 yr)−1)
No Reddening	P(AV) = 0	${79.67}_{-5.53}^{+4.51}$	$\lt {0.47}_{-0.03}^{+0.04}$	$\lt {0.93}_{-0.18}^{+0.16}$	$\lt {0.94}_{-0.37}^{+0.38}$
Ia Reddening	$P({A}_{V})={e}^{-{A}_{V}/0.334}$	${75.07}_{-6.56}^{+5.35}$	$\lt {0.50}_{-0.04}^{+0.05}$	$\lt {0.99}_{-0.15}^{+0.19}$	$\lt {1.00}_{-0.36}^{+0.43}$
High Reddening	$P({A}_{V})=2\times {e}^{-{A}_{V}/0.334}$	${69.46}_{-7.28}^{+6.15}$	$\lt {0.55}_{-0.05}^{+0.07}$	$\lt {1.09}_{-0.18}^{+0.24}$	$\lt {1.10}_{-0.40}^{+0.51}$

Notes.

^aDLT40 only detected DLT17ck because of the LIGO detection and subsequent localization; therefore, it is not considered in our rate calculations, which we report here as 95% confidence level Poisson single-sided upper limits, given zero events (Gehrels 1986). ^bWe converted the rates in units of SNuB to volumetric rates with luminosity density: (1.98 ± 0.16) × 10⁻² × 10¹⁰ ${L}_{\odot }^{B}$ Mpc³ (Kopparapu et al. 2008). ^cThe total B-band luminosity of the MW is quite uncertain; we adopt (2.0 ± 0.6) × 10¹⁰ ${L}_{\odot }^{B}$ (van der Kruit 1987).

Download table as:  ASCIITypeset image

While searching for SNe, DLT40 also reacted to LIGO/Virgo triggers during the O2 observing run, prioritizing the galaxies from the DLT40 catalog within the LIGO/Virgo localization region for each trigger. Following the LIGO/Virgo trigger of GW170817 (LIGO Scientific Collaboration & Virgo Collaboration 2017a, 2017b), DLT40 independently discovered and monitored the evolution of the kilonova AT 2017gfo/DLT17ck (Valenti et al. 2017). Given the daily cadence of the search, DLT40 is well suited to discover similar fast kilonova-like transients. In the particular case of DLT17ck, we had to revise our target priority list because the GW localization placed it near the horizon at sunset in Chile. During the course of the normal survey, however, any other kilonova-like transient would have been visible in the DLT40 search fields, given that DLT17ck itself was ∼1.5 mag brighter than our typical detection limit, out to the border of our D ∼ 40 Mpc pointed search.

One approach to measuring the rate of an astronomical transient makes use of the control time concept (Zwicky 1942; Cappellaro et al. 1993, 1997). For each ith galaxy, the control time (ct_i) is defined as the time during which a hypothetical transient is above the detection limit. It depends on the magnitude limit of each observation and the light curve of the transient. The total control time per unit luminosity of our search is computed by multiplying the ct_i control time by the luminosity of the ith galaxy and then summing over all of the galaxies in our sample⁸ :

$$\begin{eqnarray*}&&{ct}=\displaystyle \sum _{i=1}^{n}{L}_{i}^{* }{{ct}}_{i}.\end{eqnarray*}$$

The ratio between the number of transients detected and the sum of the control times for all galaxies observed gives immediately the rate as

$$\begin{eqnarray*}&&r=\displaystyle \frac{N}{{ct}}.\end{eqnarray*}$$

In order to measure the control time, the first step is evaluate the transient detection efficiency for each image or, in other words, to measure the apparent magnitude limit for transient detection. In order to do that, we performed artificial star experiments for a subset of frames, implanting stars with different magnitudes using the proper point-spread functions (PSFs) and registering the fraction of artificial stars automatically identified by our pipeline on the difference images.

Hereafter, we will adopt the magnitude corresponding to a 50% detection efficiency as the limiting magnitude for the DLT40 survey, while we use the 16% and 84% detection efficiencies as lower/upper limits to determine its uncertainty. We found that the magnitude limit of our search is in the range M_r ∼ 18–20 mag (see the left panel of Figure 2) depending on weather and seeing conditions of the specific observation. Since artificial star experiments are time consuming, instead of repeating the simulation for all of the ∼120,000 frames observed so far, we exploited a linear relation between the limiting magnitude for transient detection computed through artificial star experiments and the limiting magnitude for stellar source detection computed for each target frame (i.e., not the difference image). The latter was derived through an analytic equation using information on the noise and photometric calibration for each image. The comparisons between the two limiting magnitudes is shown in the right panel of Figure 2. In general, the limiting magnitude computed with the analytic function on the target image (y axes) is ∼1 mag deeper than the limit magnitude from artificial stars experiment (x axes). This is expected since the difference imaging technique effectively adds the template image noise to that of the target image.

Zoom In Zoom Out Reset image size

The second ingredient to measure the control time is the simulation of kilonova light curves in the time window each galaxy was observed. The time that the transient is above our detection limit contributes to the control time. The observed light curve of AT 2017gfo/DLT17ck was used as a reference, scaled to the distance of each galaxy with an explosion epoch randomly distributed in the observed time window. We took into account that kilonovae may have a range of absolute magnitudes and that they may experience a variety of host galaxy extinction due to dust. For the range in kilonova magnitudes, we varied the absolute magnitude of the kilonova using a Gaussian distribution centered on the absolute magnitude of AT 2017gfo/DLT17ck and a σ of 0.5 mag (e.g., 95% of simulated light curves have an absolute magnitude within ±1 mag of AT 2017gfo/DLT17ck). For the extinction distribution, we notice that the host environment of neutron stars mergers is often compared to the host environment of SNe Ia since both types of systems are found in early-type and star-forming galaxies (Fong et al. 2013). For this reason, we adopted for the extinction distribution $P({A}_{V})={e}^{-{A}_{V}/{\tau }_{V}}$, with τ_V = 0.334 ± 0.088 mag (Kessler et al. 2009), which we label the "SN Ia extinction" scenario. We also computed the control time using either no extinction (low extinction scenario) or an extinction distribution with a τ value inflated by a factor of 2 (high extinction scenario). We want to stress that, given that DLT17ck is the first clear case of a kilonova, any choice of absolute magnitude range and reddening law is somehow arbitrary and those quantities will be better constrained when a larger number of kilonovae are discovered.

In summary, for each galaxy, we have simulated 20,000 light curves randomly distributed in the 13 months of the search, with a range of absolute magnitudes and reddening. If at any epoch of observation the simulated light curve was brighter than our detection limit, the simulated transient would have been detected. The fraction of detected simulated transients, multiplied by the time window each galaxy was observed, gives the control time. The uncertainty on the detection limits (right panel of Figure 2) are reported as systematic errors, while the three extinction distributions used (low, similar to SNe Ia and high extinction) are reported separately.

During the 13 months of the search, the average number of observed frames per galaxy was 64, while the average control time per galaxy was 80 days. This means that any fast-evolving transients like AT 2017gfo/DLT17ck would likely not be detected a second time if the survey cadence was 2 days or larger. Our strategy of triggering a confirmation image for each new target within a few hours of the first detection well fits the need for these fast transients. Excluding AT 2017gfo/DLT17ck, which was discovered only thanks to the LIGO/Virgo trigger, no other transient with a similar fast evolution was detected. We infer a limit to the rate of kilonovae of $\lt 0.47{}_{-0.03}^{+0.04}\,\mathrm{SNuB}$ ⁹ (low extinction), $\lt 0.50{}_{-0.04}^{+0.05}\,\mathrm{SNuB}$ (SNe Ia extinction), and $\lt 0.55{}_{-0.05}^{+0.07}\,\mathrm{SNuB}$ (high extinction), where the rate has been normalized to the galaxy integrated luminosity. This translates to a limit in our Galaxy of $\lt 0.94{}_{-0.37}^{+0.38}$ (low extinction), $\lt 1.00{}_{-0.36}^{+0.43}$ (SN Ia extinction), and $\lt 1.10{}_{-0.40}^{+0.51}$ (high extinction) per 100 years. These limits and the systematic error are reported in Table 2. As a cross check, we have also computed from DLT40 the SN rates for SNe Ia, Ibc, and II that will be reported in a dedicated paper (S. Yang et al. 2017, in preparation). We stress that our SN rate estimates are consistent with previous measurements (Cappellaro et al. 1993, 1997; Leaman et al. 2011), despite the poor statistics from a few simplifications in the calculation of the control time.

In this Letter, we have used the observed light curve of a kilonova to constrain the rate of BNS mergers to less than $0.47{}_{-0.03}^{+0.04}\,\mathrm{SNuB}$ (low extinction), $0.50{}_{-0.04}^{+0.05}\,\mathrm{SNuB}$ (SNe Ia extinction), and $0.55{}_{-0.05}^{+0.07}\,\mathrm{SNuB}$ (high extinction). Since some published measurements of the BNS coalescence rates are expressed in units of comoving volume, we convert SNu rates to volumetric rates similarly to Botticella et al. (2012), that is, multiplying the SNuB rate by the galaxy B-band luminosity density reported in Kopparapu et al. (2008), $(1.98\pm 0.16\,)\times {10}^{-2}\times {10}^{10}\,{L}_{\odot }^{B}\,{\mathrm{Mpc}}^{3}$. The kilonova volumetric rate upper limit is $0.93{}_{-0.18}^{+0.16}\,{10}^{-4}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$, $0.99{}_{-0.15}^{+0.19}\,{10}^{-4}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$, or $1.09{}_{-0.18}^{+0.28}\,{10}^{-4}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (depending on the extinction law used) and is compared with previous measurements in Figure 3. Our rate is one order of magnitude higher than the BNS merger rate limit obtained by LIGO/Virgo during the observing run O1 (Abadie et al. 2010) and two orders of magnitude higher than the optimistic rates of short GRBs (Coward et al. 2012; Petrillo et al. 2013; Siellez et al. 2014).

Zoom In Zoom Out Reset image size

We can also investigate how long it would on average take for our search to discover (independently from LIGO/Virgo) a kilonova. During the LIGO O2 run (∼1 year), 117 days of simultaneous LIGO-detector observing time has been used to discover one BNS coalescence (Abbott et al. 2017a), which means there are 1/(117/365) = 3.12 BNS sources in the LIGO searching volume, while our control time for kilonovae in the same period (monitoring galaxies within 40 Mpc) is 0.22 year (on average 80 days per year per galaxy). Comparing the total luminosity of the DLT40 galaxy sample and the total luminosity of the GWGC catalog, within 40 Mpc, gives 60% of the GWGC catalog sample monitored by the DLT40 survey. In order to independently discover a kilonova we would need to run the DLT40 for 3.12/(control time * completeness) × the volume ratio of the two surveys aLIGO/aVirgo and DLT40. During the O2 run, aLIGO/aVirgo were sensitive up to a volume of 78.5 Mpc (Abbott et al. 2016b) and taking into account the different volumes of the two surveys (78.5/40)³, we would need to run DLT40 for ∼18.4 years in order to independently discover a kilonova. This explains why historical optical searches (like the Lick SN search; Leaman et al. 2011) never detected a kilonova.

Given that it is quite unlikely to independently discover a kilonova with a search like DLT40, we may at least understand what a DLT40-like search may be able to do during the O3 aLIGO/aVirgo run. During O3, LIGO will be able to detect BNS coalescences out to 150 Mpc, while Virgo should be sensitive out to 65–85 Mpc (Abbott et al. 2016a). If all kilonovae were as bright as DLT17ck, with the current DLT40 observing strategy, we could detect kilonovae within a distance of 70 Mpc. In order to cover the full Virgo volume (85 Mpc), we would need to go ∼0.4 mag deeper (to a limiting magnitude of ∼19.4 mag), hence increasing the exposure time by a factor of 2.2 (100 s per exposure, instead of the current 45 s). As DLT40 currently observes 400–600 galaxies per night with 45 s exposures, increasing the exposure time to 100 s would still allow us to observe ∼230 galaxies during a single night. Ran-domly selecting galaxies within 85 Mpc from the GLADE¹⁰ catalog in typical aLIGO/aVirgo regions (30 deg²), the average number of galaxies is ∼230—almost exactly the same number of galaxies observable by DLT40 each night with an exposure time of 100 s. Here, we neglect that the GLADE catalog is only ∼85%–90% complete in the integrated luminosity up to 85 Mpc (GLADE catalog).

Therefore, within 85 Mpc, small telescopes can still play a useful role (unless DLT17ck turns out to be a particularly bright kilonova), but the incompleteness of the available catalogs, especially for faint galaxies, may suggest that a wide-FoV strategy to directly tile the full aLIGO/aVirgo localization region may be preferred to avoid possible biases in sampling of the stellar population. In this respect, the association of GRBs (Savaglio et al. 2009) and super luminous supernovae (SLSN; Perley et al. 2016) with dwarf galaxies is a lesson learned.

Research by D.J.S. and L.T. is supported by NSF grants AST-1412504 and AST-1517649. The structure of the DLT40 web pages were developed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293 A.C. acknowledges support from the NSF award #1455090 "CAREER: Radio and gravitational-wave emission from the largest explosions since the Big Bang." The work of S.Y. was supported by the China Scholarship Council (No. 201506040044).