Optimizing Cadences with Realistic Light-curve Filtering for Serendipitous Kilonova Discovery with Vera Rubin Observatory

To assess kilonova detectability in different cadence simulations, we employed a number of metrics. We improved the existing TDEsPopMetric, designed to inject and recover diverse populations of tidal disruption event light curves, by adding the possibility to inject synthetic transients distributed uniformly in volume (with numbers increasing as a function of distance to the third power), rather than placed at a fixed distance. Light curves at a larger distance share the same properties of those at lower distances, but their apparent luminosity is fainter, making them detectable for shorter times and only when the images' magnitude limits are particularly deep.

The metrics most relevant to this work, all included as functions in the kneMetrics code, are:

• 1.  
  multi_detect: ≥2 detections >5σ;
• 2.  
  ztfrest_simple: metric that reproduces a discovery algorithm similar ¹⁸ to ZTFReST, in which sources found to be rising faster than 1 mag day⁻¹ and fading faster than 0.3 mag day⁻¹ are selected;
• 3.  
  ztfrest_simple_red: same as ztfrest_simple, but applied only to izy bands;
• 4.  
  ztfrest_simple_blue: same as ztfrest_simple, but applied only to ugr bands.

The metrics were deliberately designed to range from standard transient detection (with ≥2 detections), which typically provides only spatial information on the celestial coordinates of a source, to methods more likely to lead to source characterization—in other words, kilonova discovery. Simple detection can be crucial during gravitational-wave follow-up, but is of little use during fast-transient searches in the regular survey, especially for transients at large distances. Importantly, the conclusions of this study can be applied to a range of fast transients, including, for example, GRB afterglows and fast blue optical transients, for which light-curve sampling with spacing between 1 hr and 1 day is crucial.

There are a variety of methods in the literature to promptly identify fast-transient candidates. For example, methods are being developed for early transient classification via machine-learning techniques (e.g., Muthukrishna et al. 2019), or as part of the Photometric LSST Astronomical Time Series Classification Challenge (Kessler et al. 2019). Prior work on detecting and identifying fast transients in Rubin LSST Opsims includes Bianco et al. (2019). A simple but effective strategy to identify transients with rapidly fading or rising light curves can be based on magnitude rise- and decay-rate measurements. In this work, we consider significant fading rates to be those faster than 0.3 mag day⁻¹, which is the threshold used in real time by the ZTFReST pipeline and is expected to be particularly suitable for the discovery of kilonovae from BNS mergers (Figure 2), or rising rates faster than 1 mag day⁻¹, which can separate rapidly evolving transients from most supernovae. Within ZTFReST, this has greatly helped to separate fast-transient candidates from slower, "contaminant" sources, with ∼30% purity in "archival" data searches when considering only fade rates and thresholds tailored for each band.

Zoom In Zoom Out Reset image size

In this work, we considered kilonova models from the grid generated with the three-dimensional radiation transfer simulation code POSSIS (Bulla 2019). The model grid allowed us to explore a diversity of intrinsic properties, such as ejecta masses, as well as different viewing angles to the system.

First, we injected synthetic light curves using a single model: the GW170817-like kilonova (dynamical-ejecta mass M_dyn = 0.005M_⊙, disk-wind mass M_wind = 0.050M_⊙, and viewing angle ι = 25.8°). A half-opening angle ϕ = 30° for the lanthanide-rich region is used for this model and all the other models considered in this work. Second, we injected a population of kilonovae with the same ejecta masses of the GW170817-like model but viewed from 11 viewing angles, uniformly distributed in $\cos (\iota )$. Third, we explored kilonova detectability in an optimistic and a pessimistic scenarios, in which the kilonova properties make it particularly bright or dim, respectively. Ejecta masses were chosen to be physically realistic as determined by numerical relativity simulations, with M_dyn = 0.020M_⊙, M_wind = 0.090M_⊙ for the optimistic scenario, and M_dyn = 0.005M_⊙, M_wind = 0.010M_⊙ for the pessimistic scenario.

Cadences were made available in several releases and were grouped into "families", in which ideas that deviate from the baseline cadence were implemented and encompass parameter variations, for example in the area of the footprint. Detailed information about simulations can be found in official Rubin Observatory online resources. ¹⁹

Figure 3 shows results obtained by injecting 5 × 10⁵ synthetic GW170817-like kilonova light curves, uniformly distributed in volume out to a luminosity distance of 1.5 Gpc, in OpSim simulated cadences part of the v1.5 ²⁰ (bottom panels) and v1.7 ²¹ and v1.7.1 ²² releases (top panels). In all panels of Figure 3, each point indicates the ratio between the mean recovery fraction of an individual cadence and the maximum recovery fraction of the baseline cadence that returned the most kilonova detections. Results from those simulations are grouped by cadence family, such that all the results from cadences belonging to a certain family line up at the same abscissa. The maximum and minimum normalized recovery fractions for each family are marked by red and black triangles, respectively, and the mean recovery fraction for each cadence family is indicated by a blue circular marker to guide the eye. The results displayed in Figure 3 were obtained using the multi_detect (right panels) and the ztfrest_simple (left panels) metrics described in Section 2.1.

Zoom In Zoom Out Reset image size

In all simulations, the best baseline cadence entails individual 30 s exposures (baseline_nexp1). The baseline cadence with 2 × 15 s snaps (baseline_nexp2) performs consistently worse. The number of recovered kilonovae in cadence families simulated as part of the v1.5 cadence release (bottom panels) are relatively similar, with results comparable with the best baseline cadence within 15%. When looking at v1.7 cadences, it is evident that the best baseline performs distinctly better than any other cadence in terms of kilonova detection (multi_detect metric; top-right panel of Figure 3). The baseline cadence does a better job than most cadence families, also according to the ztfrest_simple metric. We found that rolling cadences, in which a smaller fraction of the footprint is observed in each season at higher cadence, perform significantly (∼50%–60%) worse as coded for the v1.7 release than the baseline cadence. ²³ However, rolling cadences as part of the v1.7.1 release, indicated as "new_rolling" in the figure, perform up to ∼20% better than the baseline cadence (in the figure, uncertainties are in the order of 5%–10%). In order to comparebaseline and rolling cadences with a higher statistical significance, we ran simulations in which the number of injected sources was increased to 10⁶, using a variety of surrogate models. A summary of the results can be found in Table 1.

Table 1. Kilonova Recovery Efficiencies (), Calculated with a Number of Metrics, for the Best Baseline Cadence and the Best Rolling Cadence

Kilonova	Metric	Baseline	Rolling	NKN	NKN	NKN,min	NKN,min	NKN,max	NKN,max
Model		(× 104)	(× 104)	Baseline	Rolling	Baseline	Rolling	Baseline	Rolling
GW170817	multi_detect	60.5 ± 0.8	43.1 ± 0.7	130	93	32	23	334	238
Mdyn = 0.005M⊙	blue_color_detect	6.8 ± 0.3	6.4 ± 0.2	15	14	4	3	38	36
Mwind = 0.050M⊙	red_color_detect	2.9 ± 0.2	3.1 ± 0.2	6	7	1	2	17	18
	ztfrest_simple	5.2 ± 0.2	5.6 ± 0.2	11	12	3	3	29	32
	ztfrest_simple_blue	3.8 ± 0.2	4.4 ± 0.2	8	9	2	2	22	25
	ztfrest_simple_red	1.7 ± 0.1	2.1 ± 0.1	4	4	1	1	10	12
GW170817	multi_detect	1306.3 ± 3.6	931.2 ± 3.0	40	28	10	7	101	72
<300 Mpc	blue_color_detect	141.0 ± 1.2	143.8 ± 1.2	4	4	1	1	11	11
	red_color_detect	369.5 ± 1.9	312.0 ± 1.8	11	10	3	2	29	24
	ztfrest_simple	400.3 ± 2.0	334.1 ± 1.8	12	10	3	3	31	26
	ztfrest_simple_blue	176.7 ± 1.3	173.0 ± 1.3	5	5	1	1	14	13
	ztfrest_simple_red	272.9 ± 1.6	244.5 ± 1.6	8	7	2	2	21	19
GW170817	multi_detect	31.8 ± 0.6	22.7 ± 0.5	68	49	17	12	176	127
viewing	blue_color_detect	3.4 ± 0.2	3.4 ± 0.2	7	7	2	2	20	19
angles	red_color_detect	1.4 ± 0.1	1.7 ± 0.1	3	4	1	1	8	10
	ztfrest_simple	2.6 ± 0.2	2.9 ± 0.2	6	6	1	1	15	17
	ztfrest_simple_blue	1.8 ± 0.1	2.3 ± 0.1	4	5	1	1	10	13
	ztfrest_simple_red	1.0 ± 0.1	1.1 ± 0.1	2	2	0	1	6	7
Pessimistic	multi_detect	8.9 ± 0.3	6.5 ± 0.2	19	14	5	3	50	37
Mdyn = 0.005M⊙	blue_color_detect	0.8 ± 0.1	0.9 ± 0.1	2	2	0	0	5	5
Mwind = 0.010M⊙	red_color_detect	0.3 ± 0.1	0.4 ± 0.1	1	1	0	0	2	3
	ztfrest_simple	0.5 ± 0.1	0.6 ± 0.1	1	1	0	0	3	4
	ztfrest_simple_blue	0.4 ± 0.1	0.4 ± 0.1	1	1	0	0	2	2
	ztfrest_simple_red	0.2 ± 0.1	0.3 ± 0.1	0	1	0	0	1	2
Optimistic	multi_detect	116.6 ± 1.1	87.0 ± 0.9	251	187	62	46	641	479
Mdyn = 0.020M⊙	blue_color_detect	11.9 ± 0.3	12.9 ± 0.4	26	28	6	7	67	72
Mwind = 0.090M⊙	red_color_detect	5.2 ± 0.2	5.9 ± 0.2	11	13	3	3	29	34
	ztfrest_simple	9.2 ± 0.3	10.8 ± 0.3	20	23	5	6	52	61
	ztfrest_simple_blue	6.7 ± 0.3	8.3 ± 0.3	14	18	3	4	38	47
	ztfrest_simple_red	3.2 ± 0.2	4.2 ± 0.2	7	9	2	2	18	24

Note. The efficiency was then converted into the number of expected kilonovae using the BNS merger rate $R={320}_{-240}^{+490}$ Gpc⁻³ yr⁻¹ from the GWTC-2 catalog (The LIGO Scientific Collaboration et al. 2021), where N_KN corresponds to the median rate and N_KN,min, N_KN,max correspond to the 90% symmetric credible intervals, taking the uncertainty in into account. A duration of 10 yr for the WFD survey was assumed.

Download table as:  ASCIITypeset image

When we injected GW170817-like kilonovae, the baseline cadence performed better than the best rolling cadence ²⁴ at any distance with the multi_detect metric (Figure 4, central panel), but the rolling cadence outperforms the baseline cadence beyond ∼400 Mpc with the ztfrest_simple metric. This means that a rolling cadence could enable the identification of a few more kilonova candidates than the baseline cadence at large distances. In total, 32–334 kilonovae can be expected to be detectable with the baseline cadence and 23–238 with the rolling cadence, assuming that all kilonovae are similar to the observed GW170817. However, the number of kilonovae recognizable to be fast transients (ztfrest_simple metric) would be 3–29 and 3–32 for the baseline and rolling cadences, respectively.

Zoom In Zoom Out Reset image size

Nearby kilonovae have the potential to be recognized sooner, associated with hosts at known redshifts, and can be better characterized with follow-up observations than distant (fainter) sources. To better explore detectability of nearby kilonovae, we injected 10⁶ GW170817-like kilonovae uniformly distributed in volume out to 300 Mpc, which is within the distance range where the best Wide Fast Deep (WFD) baseline cadence performs better than any of the rolling cadences simulated so far. With the best baseline cadence, we can expect up to 101 kilonovae to be detectable at at least twice in this distance range and up to 31 could be recognized to be fast-fading in at least one band. In the simulations, about 68% of kilonovae were found to be fast-fading in red izy bands (ztfrest_simple_red metric) and 44% in blue ugr bands (ztfrest_simple_blue metric). Only 37% of kilonovae found in izy bands were detected at least four times in ugr bands (Figure 4, bottom panel). The combination of transient detection, color information, and possible association with a cataloged nearby galaxy (which enables an estimate of the transient's absolute magnitude, expected to be fainter than M ∼ −16 for most kilonovae) can lead to the identification of solid kilonova candidates. For the fraction of events that are relatively nearby (below 300 Mpc), they can be followed up spectroscopically with ≳8 m class telescopes such as the Very Large Telescope, Gemini, Keck, or with the upcoming SoXS at ESO New Technology Telescope (NTT), which was designed specifically for LSST transient classification, to be classified (Schipani et al. 2016). In summary, our analysis suggests that employing more observations in redder bands is preferred to maximize scientific return.

Multimessenger observations of GW170817 constrained the viewing angle to be $\iota ={32}_{-13}^{+10}$° (Finstad et al. 2018, see also Dhawan et al. 2020). Superluminal motion from radio observations suggests a lower value for the viewing angle of about 15°–20° (Mooley et al. 2018; Ghirlanda et al. 2019).

However, merging BNS systems can be oriented in any direction with respect to the observer. We also compared the performance of baseline and rolling cadences by injecting synthetic light curves, from the grid of kilonova models simulated with POSSIS, with the same intrinsic parameters as the GW170817-like model (Section 2.2), but at different viewing angles. According to our simulations, up to 15 (17) kilonovae should be identified as fast transients in the baseline (rolling) cadence, while up to 176 (127) kilonovae should be detectable at least twice.

Finally, we assessed kilonova detectability for "pessimistic" and "optimistic" kilonova models, in which the ejecta masses make the optical emission particularly faint or bright (see Section 2.2). For the pessimistic case, only a handful of kilonovae are expected to be present in Rubin images, with at most five kilonovae expected to be recognizable as fast transients. On the other hand, the optimistic scenario could result in >50 kilonovae found to evolve rapidly with the baseline cadence and >60 with the currently best rolling cadence. A better understanding of the kilonova luminosity function is required to set more precise serendipitous kilonova discovery expectations.

Our results can be compared with a sample of related work on prospects for serendipitous detection of kilonovae by the Rubin Observatory. One difference tends to be, in some cases, the BNS merger rate assumed, for example, R = 10³ Gpc⁻³ yr⁻¹ (e.g., Scolnic et al. 2018; Setzer et al. 2019) was adopted, or up to R = 1.5 × 10³ Gpc⁻³ yr⁻¹ (Cowperthwaite et al. 2019). These rates were consistent with gravitational-wave observations before the third LIGO–Virgo–KAGRA observing run (O3), but is larger than the upper bound of the most recent rate from the GWTC-2 catalog, used in this work, of 810 Gpc⁻³ yr⁻¹ (The LIGO Scientific Collaboration et al. 2021). We also caution that the newer simulated cadences underwent some changes with respect to those available in 2018.

For a uniform population of GW170817-like kilonovae, Setzer et al. (2019) expect about 58 kilonovae to be detected during the survey using the selection criteria developed by Scolnic et al. (2018), which aim at identifying transients with multiband detections that evolve faster than most supernovae. Those results, rescaled to a lower BNS merger rate, lie well within the range of 32–334 kilonovae that could be found with the looser multi_detect metric, while they exceed the expected 3–29 kilonovae that could be found to be rapidly evolving with the stricter ztfrest_simple metric. The ∼27 expected kilonova detections from a representative population of BNS mergers in the baseline cadence, corrected for the merger rate, lies between the results that we obtained for pessimistic and optimistic kilonova models with the ztfrest_simple metric. Overall, our ztfrest_simple metric results are closer to those obtained by Cowperthwaite et al. (2019), who employed stricter selection criteria than Scolnic et al. (2018) and Setzer et al. (2019). This highlights the importance of realistic filters when estimating expected rates.