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

Current and future optical and near-infrared wide-field surveys have the potential of finding kilonovae, the optical and infrared counterparts to neutron star mergers, independently of gravitational-wave or high-energy gamma-ray burst triggers. The ability to discover fast and faint transients such as kilonovae largely depends on the area observed, the depth of those observations, the number of re-visits per field in a given time frame, and the filters adopted by the survey; it also depends on the ability to perform rapid follow-up observations to confirm the nature of the transients. In this work, we assess kilonova detectability in existing simulations of the LSST strategy for the Vera C. Rubin Wide Fast Deep survey, with focus on comparing rolling to baseline cadences. Although currently available cadences can enable the detection of more than 300 kilonovae out to 1400 Mpc over the ten-year survey, we can expect only 3-32 kilonovae similar to GW170817 to be recognizable as fast-evolving transients. We also explore the detectability of kilonovae over the plausible parameter space, focusing on viewing angle and ejecta masses. We find that observations in redder izy bands are crucial for identification of nearby (within 300 Mpc) kilonovae that could be spectroscopically classified more easily than more distant sources. Rubin's potential for serendipitous kilonova discovery could be increased by gain of efficiency with the employment of individual 30s exposures (as opposed to 2x15s snap pairs), with the addition of red-band observations coupled with same-night observations in g- or r-bands, and possibly with further development of a new rolling-cadence strategy.

Dynamical ejecta (e.g., Hotokezaka et al. 2013;Bauswein et al. 2013;Dietrich & Ujevic 2017), which arise from tidal stripping of the neutron star(s) and the neutron stars contact interface, and post-merger ejecta (e.g., Metzger et al. 2008;Fernández et al. 2015;Siegel & Metzger 2018;Fernández et al. 2019), which arise from accretion disk winds surrounding the remnant object, are characterized by low electron fractions. This scenario favors the production of heavy elements such as lanthanides and actinides via rapid neutron capture (known as the r-process), and the decay of these unstable nuclei powers the optical/infrared kilonova (e.g., Lattimer & Schramm 1974;Kasen et al. 2013; ; Barnes et al. 2016;Kasen et al. 2017).
Questions about the sources of heavy element production in the Universe and diversity in the ejecta of the kilonova population can only be answered by the detection and characterization of a large sample of sources. Unveiling such a population is difficult because kilonovae are rare (< 1 % of the core collapse supernova rate), fast (fading 0.5 mag per day in the optical), and faint transients (M −16 at peak), and hence are particularly hard to discover. Signatures of kilonovae are mostly found during the follow-up of short GRBs (e.g., Perley et al. 2009) and the follow-up of LIGO/Virgo candidates, although only for GW170817 has a counterpart been identified so far. Rates of BNS mergers are still highly uncertain, with R = 80 − 810 Gpc −3 yr −1 based on GW observations (The LIGO Scientific Collaboration et al. 2020); empirical limits on kilonovae rates by optical surveys are nearing the upper end of the gravitational-wave measurements (Andreoni et al. 2020;Andreoni & Coughlin et al. 2021).
The advent of Vera C. Rubin Observatory (Ivezić et al. 2019) presents us with a great opportunity to identify a population of kilonovae independent of any gravitational-wave or gamma-ray burst trigger, thanks to the unprecedented volume that the Legacy Survey of Space and Time (LSST) will be able to probe (see e.g., Andreoni et al. 2019;Setzer et al. 2019). Unfortunately, due to their fast fading and intrinsically underluminous properties, "detection" is not enough; it is imperative that kilonova candidates found by Rubin Observatory are recognized as such in real time so that follow-up observations can confirm their nature.
Projects exist that are dedicated to fast transient discovery in current wide-field surveys such as the Zwicky Transient Facility (ZTF; Bellm et al. 2019;Graham et al. 2019). The "ZTF Realtime Search and Triggering" ("ZTFReST" 1 , Andreoni & Coughlin et al. 2021) project, for example, employs i) alert queries, ii) forced pointspread-function photometry, and iii) nightly light curve stacking in flux space to discover fast-evolving transients such as kilonovae. ZTFReST is proving to be very effective at identifying extragalactic fast transients, having already revealed seven serendipitously-discovered GRB afterglows and at least two supernova shock breakouts in 2020 and in the first three months of 2021 (Andreoni & Coughlin et al. 2021).
In this work, we used the most recent OpSim simulations and a set of new metrics to assess the effectiveness of cadence options for un-triggered, or "serendipitous," kilonova discovery. We employed metrics that both assess Rubin Observatory's ability to simply detect the transients, as well as metrics that are designed to identify a transient as "fast" based on its flux evolution. We argue that the latter is the most appropriate metric for potentially maximizing the science output from these rare objects, and we provide suggested cadences based on this metric.

METHODS
We used the new kneMetrics 2 to recover synthetic kilonova light curves injected in OpSim simulations. The synthetic light curves are taken from Dietrich et al. (2020), which rely on the radiative transfer code POSSIS (Bulla 2019), which vary four parameters over physically viable priors: the dynamic ejecta mass M dyn , the disk wind ejecta mass M wind , the half opening angle of the lanthanide-rich dynamical-ejecta component φ and the viewing angle ι (see Dietrich et al. 2020, for more details about the adopted geometry). Examples of synthetic light curves injected in the Rubin baseline cadence can be found in Fig. 1  Observations from 30 days before peak to 30 days after peak are presented. The light curves were uniformly distributed in volume and uniformly distributed in time throughout the 10-year survey. Circles indicate the detections, solid lines show the simulated light curves in bands where at least one detection is present, and triangles indicate 5σ upper limits. The luminosity distance at which each light curve is places is also indicated. 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: • multi detect: ≥ 2 detections > 5σ • ztfrest simple: metric that reproduces a discovery algorithm similar 3 to ZTFReST, in which sources found to be rising faster than 1 mag day −1 and fading faster than 0.3 mag day −1 are selected • ztfrest simple red: same as ztfrest simple, but applied only to izy bands • 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 (FBOTS), for which light curve sampling with spacing between one hour and one 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 (PLAsTiCC; Kessler et al. 2019). Prior work on detecting and identifying fast transients in Rubin LSST Opsims include 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 −1 , 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 (Fig. 2), or rising rates faster than 1 mag day −1 , 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.

Kilonova models
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 eleven viewing angles, uniformly distributed in cos(ι). 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.

GW170817-like kilonova light curves
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 4 . Fig. 3 shows results obtained by injecting 5 × 10 5 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 5 (bottom panels) and v1.7 6 and v1.7.1 7 releases (upper panels). The results displayed in Fig. 3 were obtained using the multi detect and the ztfrest simple metrics described in §2.1. In all simulations, the best baseline cadence entails individual 30 s exposures (baseline nexp1). The baseline cadence where 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 Fig. 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 4 for example https://github.com/lsst-pst/survey strategy 5 https://github.com/lsst-sims/sims featureScheduler runs1.5 6 https://github.com/lsst-sims/sims featureScheduler runs1.7 7 https://community.lsst.org/t/survey-simulations-v1-7-1-release-april-2021/ 4865 cadence, perform significantly (∼ 50 − 60%) worse as coded for the v1.7 release than the baseline cadence 8 . However, rolling cadences 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 compare baseline and rolling cadences with a higher statistical significance, we ran simulations in which the number of injected sources was increased to 10 6 , using a variety of surrogate models. A summary of the results can be found in Tab. 1. When we injected GW170817-like kilonovae, the baseline cadence performed better than the best rolling cadence 9 at any distance with the multi detect metric (Fig. 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 8 Significant changes to the Opsim approach to simulate rolling cadence strategies were implemented for v1.7, such that these simulations should be considered more reliable (see Bianco   (ztfrest simple metric) would be 3-29 and 3-32 for the baseline and rolling cadences, respectively. Nearby kilonovae have the potential of being 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 6 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 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 4 times in ugr bands (Fig. 4, bottom panel). The combination of transient detection, color informa-tion, and possible association with a catalogued 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 (VLT), 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.

Exploring the kilonova light curve parameter space
Multi-messenger observations of GW170817 constrained the viewing angle to be ι = 32 +10 −13 deg (Finstad et al. 2018), see also Dhawan et al. (2020). Superluminal motion from radio observations suggests a lower  , and a ZTFReST-like metric applied only to red (izy) bands (red lines). One million sources were injected uniformly distributed in volume between 10 Mpc and 1.5 Gpc. Center: Efficiency as a function of luminosity distance; 5 × 10 5 sources were injected at regular intervals. For the ztfrest simple metric, due to the rapidly growing rate at the edge of the sensitive volume, small differences between detection efficiencies for the rolling and the baseline cadence beyond ∼ 400 Mpc is enough to yield an improvement of ∼ 20% in kilonova detection. Bottom: Ratio between the detection efficiency in redder izy bands and in bluer ugr bands for multi-detection and ZTFReST-like metrics. Employing redder filters presents clear advantages at lower distances, where spectroscopic and multi-wavelength follow-up observations are possible.
value for the viewing angle of about 15-20 deg (Mooley et al. 2018;Ghirlanda et al. 2019). However, merging BNS systems can be oriented in any direction with respect to the observer. We compared the performance of baseline and rolling cadences also by in-jecting synthetic light curves, from the grid of kilonova models simulated with POSSIS, with the same intrinsic parameters as the GW170817-like model ( §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 §2.2). For the pessimistic case, only a handful of kilonovae are expected to be present in Rubin images, with at most 5 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 serendipidous kilonova discovery expectations.

CONCLUSION
Rubin Observatory has a great potential of revealing a population of kilonovae during the WFD survey, in addition to discoveries made following up GW triggers. We injected synthetic kilonova light curves into simulated Rubin observations to assess which ones of the available cadences can maximize serendipitous kilonova discovery. We demonstrated that, for the WFD survey, the simulated baseline cadence with single 30 s exposures should be greatly preferred over 2 × 15s consecutive snaps for kilonova discovery.
Rolling cadences are expected to be particularly suitable for fast transient discovery (e.g., Andreoni et al. 2019). We found that the development of rolling cadences has significantly improved from OpSim version v1.7 to v1.7.1. While this indicates progress, the baseline plan may still be preferred over any other cadence family currently available due to a larger efficiency at detecting more nearby (and therefore brighter) fast transients, easier to follow-up and classify with other telescopes. We recommend simulating new rolling cadences further optimizing the algorithms used in v1.7.1, possibly maximizing the exposure time in each band (barring u-band) rather than using snap pairs.
We found strong evidence that red izy bands are preferred for kilonova discovery at distances below 300 Mpc, in agreement with the results of other studies such as, for example, Almualla et al. (2021) and Sagués Carracedo et al. (2021). This is expected because kilonovae appear as red and rapidly-reddening transients due to heavy r-process elements synthesised in neutron-rich Table 1. Kilonova recovery efficiencies ( ), calculated with a number of metrics, for the best baseline cadence and the best rolling cadence. The efficiency was then converted into the number of expected kilonovae using the BNS merger rate R = 320 +490 −240 Gpc −3 yr −1 from the GWTC-2 catalog (The LIGO Scientific Collaboration et al. 2020), where NKN corresponds to the median rate and NKN,min, NKN,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.

Kilonova
Metric baseline rolling ejecta. At redder wavelengths, kilonova light curves can be brighter and longer-lived, especially if the system is viewed from equatorial viewing angles (e.g., Bulla 2019). Very rapid "blue" kilonovae could be found at larger distances ( Fig. 4) due to the greater sensitivity of g and r filters, however, these kilonovae might be rarer and more difficult to classify spectroscopically. Therefore, we recommend that the number of izy observations is increased in the WFD cadence plan. Such red-band observations would be particularly effective, scientifically, if coupled with at least one observation in g-(preferred) or r-band on the same night, so that kilonovae can be separated photometrically from other transients and their temperature evolution can be measured. A recommendation for same-night multi-band photometry in LSST has already been put forward for example by Andreoni et al. (2019) and Bianco et al. (2019). In particular, Bianco et al. (2019) address the advantages of acquiring sets of three exposures per field in the same night, in two filters and appropriately spaced in time, towards rapid identification of rare fast transients. Major uncertainty in the results of this work results from our limited understanding of the BNS merger rate and the kilonova luminosity function. Systematic kilonova searches during gravitational-wave follow-up (e.g., Kasliwal et al. 2020), short GRB follow-up (e.g., Gompertz et al. 2018;Rossi et al. 2020), and un-triggered wide-field surveys (e.g., Doctor et al. 2017;Andreoni et al. 2020;Andreoni & Coughlin et al. 2021;McBrien et al. 2021) are expected to improve those measurements significantly before Rubin Observatory's first light.