Photometric prioritization of neutron star merger candidates

Rapid identification of the optical counterparts of Neutron Star (NS) merger events discovered by gravitational wave detectors may require observing a large error region and sifting through a large number of transients to identify the object of interest. Given the expense of spectroscopic observations, a question arises: How can we utilize photometric observations for candidate prioritization, and what kinds of photometric observations are needed to achieve this goal? NS merger kilonova exhibits low ejecta mass (~5x10^-2 solar mass) and a rapidly evolving photospheric radius (with a velocity ~0.2c). As a consequence, these sources display rapid optical-flux evolution. Indeed, selection based on fast flux variations is commonly used for young supernovae and NS mergers. In this study, we leverage the best currently available flux-limited transient survey - the Zwicky Transient Facility Bright Transient Survey - to extend and quantify this approach. We focus on selecting transients detected in a 3-day cadence survey and observed at a one-day cadence. We explore their distribution in the phase space defined by g-r, g-dot, and r-dot. Our analysis demonstrates that for a significant portion of the time during the first week, the kilonova AT 2017gfo stands out in this phase space. It is important to note that this investigation is subject to various biases and challenges; nevertheless, it suggests that certain photometric observations can be leveraged to identify transients with the highest probability of being fast-evolving events. We also find that a large fraction (~0.75) of the transient candidates with |g-dot|>0.7 mag/day, are cataclysmic variables or active galactic nuclei with radio counterparts.


INTRODUCTION
The classification of astronomical transients usually requires spectroscopic resources.To date, only about 10% of the transients reported worldwide to the Transient Name Server (TNS 1 ) have spectroscopic observations (Kulkarni 2020).The alternative of photometric classification has partial success (Poznanski et al. 2002).
A related topic is how to prioritize follow-up resources for candidates of optical (ultraviolet to infrared) emission that originate from compact-object mergers and additional types of fast-evolving transients.With the increase in the gravitational-wave detector horizon, deeper searches are re-quired, which in turn will increase the number of transient candidates that are found in the gravitational-wave error region.Furthermore, in order to study the kilonova spectral evolution from early times after the merger, very rapid classification or at least target prioritization is required.
A feature of NS-merger-like transients is that they involve ejecta with: low-mass (Mej ∼ 5 × 10 −2 M⊙), opacity of < ∼ 1 cm 2 g −1 , and high velocities (vej ∼ 0.2c), where c is the speed of light (e.g.Rosswog et al. 2018;Waxman et al. 2018;Nakar 2020).In the first few days after the merger, the ejecta are optically thick, resulting in a very fast evolution of the light curve (compared to supernovae and most other transients).For example, in the case of AT 2017gfo, the optical counterpart of GW 170817 (Abbott et al. 2017), the ejecta velocity can be estimated from the photometric observations alone (see e.g., Waxman et al. 2018).Specifi-cally, when enough observations are available, and the distance is known, it is possible to fit the photospheric radius as a function of time and obtain its velocity.Here, we examine a simpler approach of using the color and flux derivatives of the transients.Using flux derivatives is not new.For example, Khazov et al. (2016) identified young core-collapse SNe based on their rise-time evolution, while Andreoni et al. (2021) applied this for NS-merger candidates.Bianco et al. (2019) solved the reverse problem and used this to estimate the minimal cadence requirements for an LSST-like survey to identify fast transients like compact-object merger afterglows, and rapidly evolving blue transients (e.g., Drout et al. 2014;Ho et al. 2021;Ho et al. 2019a;Ho et al. 2019b;Ofek et al. 2010;Ofek et al. 2021).The Bianco et al. (2019) estimate is based on simulated Type-Ia SNe, and a few observed corecollapse SNe.Here, we inspect the position of various transients, including AT 2017gfo, in the color vs. magnitude timederivatives phase space.We find that a simple combination of the transient color and the magnitude-derivatives in one or two bands is sufficient, in most cases, for removing a large fraction of unrelated transient sources and efficiently prioritizing follow-up observations.
In §2, we present the sample of background transients that we use, while in §3, we calculate the color and magnitudederivatives of objects in this sample, as well as for AT 2017gfo.In §4 we present the distribution of transients in the g − r, ġ, and ṙ phase space, in §5 we discuss the nature of the common outliers in this phase space, and we conclude in §6.Throughout the paper, the dot symbol above the band name denotes a time derivative.Throughout the paper, we use tools from Ofek (2014).

TRANSIENTS SAMPLE
To estimate the background of transients we use a sample of sources, which are part of the Zwicky Transient Facility (ZTF; Bellm et al. 2019a, Graham et al. 2019, Dekany et al. 2020, Bellm et al. 2019b) Bright Transient Survey (BTS;Fremling et al. 2020;Perley et al. 2020) found during 2018-20212 .Here, we summarize how these transients are identified in near real time.
The ZTF pipeline (Masci et al. 2019) does the image calibration and produces difference images based on the Zackay et al. (2016) algorithm.PSF-fit photometry runs on the difference images and sources in different epochs are matched to generate light curves.The pipeline generates about 10 6 alerts per night (Patterson et al. 2019, Masci et al. 2019).These are then filtered down to a small number of candidates (∼ 500 per night in 2018 (Fremling et al. 2020) and < 50 per night after that (Perley et al. 2020)).At the end of the night an astronomer reviews this list, judges which of the candidates are likely genuine astrophysical transients, and triggers spectroscopy based on the expected luminosity at peak (Fremling et al. 2020).Candidates with prior variability are discarded, such that a large fraction of active galactic nuclei (AGNs) and cataclysmic variables (CVs) are rejected.However, some are saved as transients nevertheless.Transient candidates are usually classified with the SEDMachine on the P60 telescope (Ben-Ami et al. 2012;Blagorodnova et al. 2018).
We downloaded alert photometry through the Fritz broker (Duev et al. 2019;Kasliwal et al. 2019).However, nightly alert packages are also publicly available from the ZTF Alert Archive3 for nights since June 2018.The ZTF Avro Alert package4 provides the tools to read the alerts and filter them by object ID to extract only the ones that were selected for the BTS5 sample.Alternatively, alert light curves can be obtained from different community brokers.
The BTS is based on the public three-day cadence survey, and only these observations are used to identify transients.However, many of these transients have ZTF observations obtained at a higher cadence.Here we use all the available ZTF data (i.e.not only the BTS survey data).As discussed in the next section, we are using only observations which have one day (or faster) cadence.This may introduce biases when trying to estimate the probability distribution function of fastevolving transients in the color and magnitude-derivatives phase space.For example, fast transients have a higher probability of being missed.Nevertheless, given the large number of transients in the BTS sample, including some fast-evolving transients, it is likely good enough for our objective of target prioritization and reducing the follow-up load.
BTS aims to obtain spectra for all bright transients and is 97% complete for SNe brighter than 18th magnitude, 93% for 18.5 and 75% for 19 (Perley et al. 2020).
The fact that the BTS sample magnitude limit is about two magnitudes brighter compared to the ZTF limiting magnitude is important to our purpose.Specifically, without this magnitude limit, our sample will be dominated by transients detected near maximum light, and will be biased towards objects with slow magnitude evolution (i.e., the magnitude derivative near the peak is minimal).
We corrected the light curves for the Galactic extinction (Schlegel et al. 1998), assuming RV = 3.08 (Cardelli et al. 1989).We only used photometric measurements that have errors smaller than 0.3 mag.We note that using different measurement errors cut of 0.05 mag or 0.15 mag do not changes the results significantly.We also removed sources brighter than magnitude 14.5.The final sample we use is an order of magnitude smaller than the full BTS sample.The main reason for this is our requirement for selecting objects that were observed with a 1-day cadence and small photometric errors.
The color and magnitude-derivative distribution of sources on the sky depends on when the last search of the sky was conducted.For example, a survey with one month cadence will mostly detect old transients that evolve slowly, while the new transients found by a one-day cadence survey will mostly be young transients or transients near their maximum light.For that reason, we present results for two samples.The first sample, is for all the observations, regardless of the age of the transient, and the second sample, is for only the first epoch of each unique transient.In Table 1 we list the number of magnitude-derivative measurements and unique objects for each transient class.

CALCULATION OF THE COLOR AND MAGNITUDE-TIME DERIVATIVES
For each transient in our sample, we calculate its g-and rband magnitude time derivatives ( ġ, and ṙ, respectively) and its g − r color at various epochs.In order to calculate the derivatives and colors we select light curves with at least three observations in two consecutive nights in the g and r band, respectively.For each filter, at least one of the three observations of each filter is in one night, and at least two others are in the previous or next night.Practically, this was done by searching for time windows of at least 18 hours and less than 30 hours (i.e., two consecutive nights) in which there are at least three g-band observations, and independently also, at least three r-band observations that were taken on two successive nights.The ZTF data is almost always collected with at least two observations per night.Therefore, our requirement for three data points within the time window is not diluting the sample size significantly.We selected only photometric data points with an uncertainty smaller than 0.3 magnitude.For each set of data points fulfilling these criteria, we fitted a first-degree polynomial as a function of time (Press et al. 2002;Gould 2003), separately to the g and r band data.The first-degree polynomial was fitted after subtracting the window mid-time from the time of the observations6 .The slope of the polynomial gives us ġ or ṙ, while the intersection gives us the mean magnitude, from which the color can be determined7 .We also rescaled the ZTF r-band errors by a factor of ∼ = 1.17 ∼ = (1.38).The reason for this is that when doing the linear fit for each time window, without the scaling, the average, over all the fits, of the χ 2 /dof , was about 0.96 and 1.38, for the g, and r-band, respectively.Finally, these magnitudes are corrected for Galactic extinction.We note that our ġ were measured over one day time scale, and that for fast-evolving transients the magnitude derivative over shorter time scales may be larger.Table 2 provides all the derived values and their estimated errors.In total, we have 3225 epochs for 674 unique transients.The transient types include: active galactic nuclei (AGN), cataclysmic variables (CV), luminous blue variables (LBV), superluminous SN (SLSN), Type II, IIn, Ia, Iax, Ia-91T, Ia-CSM, Ib, Ic, Ic broad line, and tidal disruption event (TDE) candidates (see Table 1).The classifications are derived from the original ZTF-BTS classifications (Fremling et al. 2020;Perley et al. 2020), which were made through human inspection of the spectra and are therefore subjective.

DISTRIBUTION OF THE COLOR AND MAGNITUDE-TIME DERIVATIVES
Figure 1 shows the g − r vs. ġ values of all BTS transients that pass our cuts, as well as the position as a function of time of AT 2017gfo.Figure 2  There are some differences between the two photometric compilations of AT 2017gfo, especially at early times < ∼ 1 day.This is presumably because it is difficult to estimate the magnitude-time-derivative when the magnitude differences between adjacent observations are small.The Waxman et al. ( 2018) compilation seems to result in smoother estimates of g − r, ġ, and ṙ.Therefore, for our objective, they are likely preferred.Regardless of the actual magnitude-derivative of AT 2017gfo, we expect the ġ and ṙ to be negative and have larger values than other transients at some early time after the merger ( < ∼ 0.5 day).To summarize, AT 2017gfo seems, as expected, to evolve faster than other transients, for the majority of the time.The two compilations demonstrate that accurate estimates of ġ, ṙ, and g − r are not straightforward to obtain.Nevertheless, even with the current uncertainties, it is possible to separate compact-object mergers from the majority of transients.
Figure 3 presents the ġ distribution for several classes of objects, including all measurements (heavy black), first measurement for each unique object (heavy blue), and selected classes of transients.The first important feature, seen in this plot, is that the first observation of each transient tends to have a more extreme ġ value.Second is that some classes of objects, most noticeably CVs, tend to produce a higher fraction of extreme ġ values.
To visually demonstrate that it is possible to separate AT 2017gfo-like events from more common transients, we would like to find a transformation of the three observables g − r, ṙ and ġ that separates the two populations well and along the entire time evolution of AT 2017gfo.For simplicity, we use a linear transformation and perform a principalcomponents analysis on the g−r, ṙ and ġ of all the AT 2017gfo observations from the Waxman et al. (2018) compilation.We find that 99.8% of the variance content is in the first and second principal components.Figure 4 shows the first principal component (PC1) vs. the second principal component (PC2) for all the transients in the BTS sample as well as AT 2017gfo.The equations for the principal components, after the extinction correction, are: P C3 = 0.6634(g − r) + 0.2659 ġ + 0.6995 ṙ.
(3) This is not necessarily the best linear transformation that separates the two populations, but it seems good enough for presentation purposes.Figure 4 demonstrates that the infor-  mation content in the g − r color is not negligible.Specifically, it is possible to separate between AT 2017gfo and the BTS transient population up to day ∼ 10, after the merger.
There are multiple ways to use this information for target selection.This includes: preferring targets that are located at low-density regions of the phase space we explore in this work, and using Table 2 to estimate the rough observed probability distribution of events in this phase space9

OUTLIERS IN Ġ
It is worthwhile to inspect the nature of the outliers seen in Figure 1.In Table 4 we list all the 47 outliers with | ġ| > 0.7 mag day −1 .An interesting finding is that about 50% of the entries in this table are CVs, and about 25% are AGNs.
In most cases, these can be easily identified by the fact that they likely had a previous flare (or large derivative phase) in the data (i.e., ∆Tg larger than a few days).The 13 events associated with AGNs are due to six unique objects, while among 24 events associated with CVs, there are 15 unique objects.
-1 -0.5 0 0.5 1 The fast rise time of dwarf novae (that likely dominates our sample) is well known.For example, Cannizzo & Mattei (1998) found that the median magnitude derivative of SS Cyg outbursts is about 2 mag day −1 , with tail extending to about 4 mag day −1 .
We also added to Table 4 a column indicating the total flux of the brightest VLASS (Lacy et al. 2020) radio source found within 3 ′′ of the transient position 10 .Lyke et al. (2020) found that about 3% of the quasars targeted by the SDSS have a radio counterpart in the FIRST catalog (e.g., Becker et al. 1995).This is in sharp contrast to the finding that all our AGN labeled events with | ġ| > 0.7 mag day −1 have a radio counterpart, brighter than the FIRST detection limit.

DISCUSSION
We analyze the distribution of transients, identified by ZTF-BTS, in the g −r, ġ, and ṙ phase space.These transients have daily observations but were selected through a 3-day cadence survey.The primary objective of this analysis is to establish criteria for prioritizing follow-up observations of GW events.
Within this sample, we identify three key challenges: (i) The sample's foundation rests on a particular observation strategy and selection criteria that might not be applicable to other sky survey approaches.(ii) The magnitude derivatives are computed on a daily timescale, whereas NS merger events (and other rapidly evolving transients) benefit from measurements taken within a span of a few hours (similar to our comparison with AT,2017gfo).(iii) Our study is built upon specific filters.
With respect to the survey cadence, we can discuss three representative scenarios for a hypothetical survey that is being used prior to the GW trigger.The first scenario is that the last epoch was 3 day prior to the GW trigger (e.g., LSSTlike cadence).In this case, the distribution of transients found in such a survey will be narrower than the distribution of the orange points in figures 1, 2, 4. The second scenario is a sky survey with a nightly cadence11 .In this case, the distribution of objects in the g − r, ġ, ṙ phase space may be wider than the distribution of the orange points in figures 1, 2, 4. For example, inspecting the ġ distribution in Figure 3, in the worst-case scenario the density of the wings of the distribution may be up to 3 times higher compared to their current location in this plot12 .This is still about two orders of magnitude below the peak of the objects in Figure 3. Therefore, in this case, transients with ġ 0.5 mag day −1 are still rare.The third scenario is a sky survey with multiple observations per night13 .In this case, we are in an unexplored territory of the transient-durations phase space.However, current limits suggest that the rate of fast transients is likely lower than the rate of transients with a few days time scale (e.g., Perley et al. 2020;Ho et al. 2023).
With respect to the second problem (i.e., magnitude derivatives based on a one-day cadence survey).Here, the problem is that we wish to make a photometric decision in a few hours of observations, rather than one day.In this case, the question is: Are the one-day-based derivative relevant for measurements on a shorter time scale?The answer is that we do not know.However, in this case, we can argue that at least14 2/3 of random transients that are found, by a 3-day cadence survey, will be older than 1 day.Therefore, in this case, one can use the distributions of objects in Figure 1,2,4, as an order of magnitude indicator.A related issue is how well can we measure magnitude derivatives on a shorter time scale.This of course depends on how rapid the variations are, and on the photometric precision.Currently, groundbased sky surveys are limited to about 1% absolute photometric accuracy (e.g., Padmanabhan et al. 2008;Ofek et al. 2012;Schlafly et al. 2012).Reaching one percent photometric accuracy for a 21-st magnitude target, using ∼ 0.5 m class telescopes is an order of magnitude more time-efficient than obtaining spectroscopy for a similar target using a 3-m class telescope.Assuming one percent accuracy per measurement, for a target with ġ > 0.5 mag day −1 , we need a time span of 2 hr.
Another limitation is the fact that our work is based on specific filters (g and r).In the first few days after the explosion, even transients with low-ejecta mass (∼ 10 −2 M⊙), are optically thick, and their spectrum is roughly described by a black-body curve.Therefore, the ġ and ṙ can be roughly interpolated and extrapolated, but the accuracy of this process is hard to quantify.A related problem is that filters (even if they have the same name), on different telescopes are not identical, which will introduce a color term.In principle, these color terms can be estimated in advance.
Although obtaining observations in multiple filters is not expensive, relative to spectroscopy, sometimes such observa-tions will not be available.In this case, we note that observations with a single filter, or even color information, without temporal information can be used, with lower effectiveness, for target prioritization.
Figure 1 demonstrates that there may be a few fastevolving transients whose ġ is significantly different from zero.However, these objects are relatively rare, and therefore, we will still achieve our goal of reducing the number of targets for follow-up.Among the most frequent fast-evolving transients are SNe found at very early times after the explosion (e.g., the green star in Figure 1; Gal-Yam et al. 2014, Yaron et al. 2017), fast transients (e.g., Drout et al. 2014, Ofek et al. 2010, Ho et al. 2021, Ofek et al. 2021), and likely also GRB afterglows (e.g., Cenko et al. 2013).The expected rate of SNe brighter than magnitude15 21 is about two orders of magnitude higher than the GRB rate.This can be estimated by multiplying the volume to which we can detect a −18 absolute magnitude SN (the median abs.mag. of SNe; Perley et al. 2020), by the SNe rate (∼ 10 −4 Mpc −3 yr −1 ; Leaman et al. 2011), and comparing it with the GRB rate multiplied by the fraction of GRBs that are visible (∼ 500 yr −1 ; Cenko et al. 2009).Furthermore, relativistic transients with no GRB emission are likely not as common as GRBs (Ho et al. 2022).
Finally, an interesting finding is that about 75% of the ġ 0.7 mag day −1 events in our sample are either AGNs or CVs.Many of these objects have recurring outbursts and can be identified in variability studies (e.g., Wozniak et al. 2002;Drake et al. 2014;Ofek et al. 2020), and be rejected as kilonova candidates.

Figure 1 .
Figure 1.The g − r color vs. the time-derivative of the g-magnitude for the BTS transient sample and AT 2017gfo, the optical afterglow of GW 170817.The thick gray line shows the observed time evolution of AT 2017gfo (Waxman et al. 2018), where the color of the filled circles, as well as the numbers, indicates the kilonova age since the NS merger (see colorbar).The thin gray line and the empty circles are the same but for the Arcavi et al. 2017 compilation.The black dots show transients from the BTS sample (Fremling et al. 2020, Perley et al. 2020).The orange dots are for the first epoch of each unique transient.The green square shows AT 2018lqh during the decay, immediately after maximum light (Ofek et al. 2021); the red box represents the M85-OT 1 luminous red nova transient during the plateau phase (Kulkarni et al. 2007, Ofek et al. 2008); the green triangle shows SN 2019hgp (Gal-Yam et al. 2022); and the green star shows the early observations of SN 2013fs taken during the light curve rise(Yaron et al. 2017).The black arrow shows the direction of the reddening with E B−V = 0.2 mag(Cardelli et al. 1989).All magnitudes were corrected for Galactic extinction.

Figure 3 .Figure 4 .
Figure 3. ġ distribution (histograms) for selected classes of transients.The black line represents the entire sample and the blue line is the first epoch of all unique transients.The histograms are normalized such their sum is unity.The heavy dashed lines shows the ġ distribution calculated using a magnitude-error cutoff of 0.1 mag instead of 0.3 mag.

Table 1 .
Number of BTS transients and photometric measurement by class.

Table 2 .
BTS transients, magnitudes and magnitude derivatives after Galactic extinction correction.δ symbols indicate uncertainties The first five lines are displayed.MJD refers to the mid-time of the g-band observations.∆Tg is the time since the first g-band detection.For repeating events this refers to the time since the first detected event.The full table is available in the electronic version of this manuscript.

Table 3 .
AT 2017gfo colors and magnitude-derivatives as a function of time.

Table 4 .
Subset of Table2, containing outliers with | ġ| > 0.7 mag day −1 .∆Tg for repeating events this refers to the time since the first detected event.Also added is the total S-band radio flux of the VLASS radio counterparts within 3 ′′ .