Data-driven Detection of Multimessenger Transients

A diffuse $\mathrm{TeV}$–$\mathrm{PeV}$ flux of astrophysical neutrinos has been discovered by IceCube (Aartsen et al. 2013). While the exact nature of the emission remains elusive, its apparent isotropy suggests that it originates from relatively weak extragalactic sources. Possible sources of ultra high-energy cosmic rays (UHECRs) and neutrinos include long- and short-duration gamma-ray bursts (GRBs), as well as core-collapse supernovae (CC-SNe) with choked jets or shock breakouts (Kashiyama et al. 2013; Senno et al. 2018). The energy density of the astrophysical neutrinos is comparable to that of the isotropic $\gamma$-ray background and to that of UHECRs (Meszaros et al. 2019). This indicates that multimessenger interpretation may lead to breakthroughs in our understanding of cosmic ray (CR) accelerators. It may elucidate the connection between supernovae (SNe) and GRBs, and may shed light on the nature of their central engine. Searches commonly take the form of either spatio-temporal clustering of neutrinos, or their direct association with steady or transient sources (Aartsen et al. 2015).

We search for clustering in two all-sky samples. The first comprises IceCube event lists of track-like muon neutrino candidates, taken between 2011 and 2012 (MJD 55694–56415; Aartsen et al. 2017). The second is of ANTARES muon neutrinos, observed within the same time period (Adrián-Martínez et al. 2014). The nominal inputs to our RNN are based on neutrino event density metrics (see Appendix A.2). The densities are defined with respect to a given location on the sky. For the IceCube sample, the data are split into four logarithmically spaced bins in the energy proxy, between 10 $\mathrm{GeV}$ and 8 $\mathrm{PeV}$. For ANTARES the metrics are inclusive in energy. The data are integrated over 24 hr time periods, which effectively avoids dependence of event rates on R.A. (Aartsen et al. 2015). The response of the IceCube and ANTARES detectors depends on the zenith of the observed event, which is therefore added as input to the RNN. In cases where an IceCube source is not within the field of view of ANTARES, we use background data instead. In total, we have $\eta =6$ inputs per time step. The collection of inputs is derived for ${\tau }_{\mathrm{RNN}}=15$ days periods. The first ${\tau }_{\mathrm{enc}}=10$ days are assumed to only contain background events. Transient signals are searched for within the next ${\tau }_{\mathrm{dec}}=5$ days interval. In total, the RNN receives $6\times 15=90$ inputs.

The RNN is trained to predict the neutrino event densities in each of the five days being probed, for a particular sky position. As part of anomaly detection, these predictions are meant to correspond to the background. Potential transient signals must therefore be removed from the training data set. This is done by scrambling neutrino events in R.A. and time of detection.

The test statistics for detection are based on differences between the predictions of the network and the true data. For our particular definitions of observables, one can not assume a priori that a particular statistical model (e.g., Poissonian counts) would be valid. We therefore calculate detection probabilities in a more general approach, using simulations. We consider the zenith to be an auxiliary RNN parameter, as it is not directly used to derive detection probabilities. Rather, we bin the data in zenith into 90 intervals of $2^\circ$ width, in which the detector response is approximately uniform. We then generate $\sim {10}^{6}$ scrambled background sequences for each bin and evaluate them with the RNN. The derived distributions of test statistics are used to parameterize the correspondence between TS- and $p$-values, accounting for all spatial and temporal trials.

We proceed to use the trained network to evaluate the original (non-scrambled) data, carrying out an all-sky grid search with 0.3° spacing. The corresponding spatial distribution of $p$-values is shown in Figure 2. The most significant position is $\{{\rm{R}}.{\rm{A}}.,\,{\rm{D}}{\rm{e}}{\rm{c}}{\rm{l}}.\}=\{{163}^{\circ },-26.5{}^{\circ }\}$. It has ${p}_{\nu }=1.9\times {10}^{-7}$ (pre-trials) and ${p}_{\nu }=0.1$ (post-trials), corresponding to $5.2\sigma$ and $1.6\sigma$ significance. The result does not correspond to any astronomical object of interest, and is consistent with a statistical fluctuation. We also use the outputs of the RNN to perform a correlation analysis between the IceCube and ANTARES events (see Appendix A.2), finding no significant result.

Zoom In Zoom Out Reset image size

Our conclusions on the existence of neutrino sources are consistent with previous studies of these data, which did not detect any source (${p}_{\nu }=0.6$, by Aartsen et al. 2015). In this case, however, the analysis is done without the need to explicitly define likelihood functions for the background or sources. The study is performed on a combined IceCube and ANTARES data set. Another advantage of our approach, is that there is no need to model the relative response between the two neutrino observatories. This example also illustrates the flexibility of the methodology regarding observables. Our choice of inputs and test statistics is motivated by the properties of the data sets, but is by no means unique. However, our framework is designed to provide self-consistent detection probabilities in general, given a set of primary (e.g., neutrino densities) and auxiliary (e.g., zenith) parameters. Finally, we note that the limitations of the public data sets constrain us to $\geqslant$ 1 day time bins. Provided that the full data-streams of the experiments become available, the same approach would be applicable for real-time searches on shorter timescales.

An alternative to auto-correlation analyses is to search for cross-correlation between different messengers. An important physical example is that of CC-SNe with relativistic outflows (Murase & Ioka 2013; Cano et al. 2017). CC-SNe may accelerate UHECRs, which produce $\gamma$-rays and neutrinos, for instance, via ${pp}$ or $p\gamma$ interactions. A direct connection between neutrinos and their SN counterparts can be made for events that also exhibit high-energy emission, e.g., GRBs. However, $\gamma$-rays are not always observed in association with CC-SNe. For example, depending on the environment within and around a source, $\gamma$-rays may become attenuated due to $\gamma \gamma$ interactions (Boncioli et al. 2019). A few percent of CC-SNe are estimated to harbor undetectable jets, which are not powerful enough to punch through their progenitors and winds. Such events are associated with choked jets, and possibly shock breakouts (Kashiyama et al. 2013). Conversely, those jets that are successful are mostly launched off-axis with respect to the observer, and thus are also undetected (Denton & Tamborra 2018). Finally, even when the high-energy emission is beamed toward Earth, individual SNe are unlikely to be identified as sources of neutrinos, as illustrated in Figures 3(A)–(B). This motivates performing a stacking analysis, combining observations from many events.

Zoom In Zoom Out Reset image size

We simulate a neutrino stacking analysis, intended to identify an accumulated neutrino over-density in spatio-temporal coincidence with CC-SNe (Senno et al. 2018). (This should be distinguished from the more direct cross-correlation approach, which involves optical follow-up of specific neutrino events; Morgan et al. 2019.) The purpose of the RNN is to provide a background model for the joint neutrino and optical observation. Optical transients are simulated according to projected observations with the Wide-Fast-Deep survey of the upcoming Vera C. Rubin Observatory (previously referred to as the Large Synoptic Survey Telescope, LSST; LSST Science Collaboration et al. 2009). These SNe are only used to define time-intervals, during which neutrino flares may occur. Neutrino signals, based on our IceCube sample, are integrated over these intervals.

The inputs to the RNN are of two types. For a given time step, the first consists of the IceCube neutrino densities (as for the previous example). The second type consists of LSST signal-to-noise metrics (Zackay & Ofek 2017) in five optical bands (see Appendix A.3). The zenith of observation is also added as an auxiliary parameter. We thus have $\eta =10$ inputs per time step in total. SNe light curves evolve over days–weeks. We therefore construct ${\tau }_{\mathrm{RNN}}=20$ days data sequences, with the first ${\tau }_{\mathrm{enc}}=10$ days representing the background, and the next ${\tau }_{\mathrm{dec}}=10$ the search period. In total the RNN receives $10\times 20=200$ inputs.

We begin by creating background samples. These represent periods of joint neutrino and optical observations, in which no transient exists in either messenger. The neutrino data are derived by the scrambling procedure described above. The optical data are chosen based on the true peak time of simulated SNe. The background samples are used to train the RNN. For a particular sky position, we then derive individual test statistics for source detection with each messenger, which are calibrated to $p$-values as a function of zenith.

In the next step of the analysis, we simulate stacked neutrino signal samples. For the source model, we assume that optically detected SNe are in fact GRBs, and that neutrinos are emitted during the short prompt phase of explosions. The fluence of muon neutrinos is modeled as ${{ \mathcal F }}_{{\nu }_{\mu }}\propto {{ \mathcal E }}_{\mathrm{CR}}/{D}_{{\rm{L}}}^{2}$, where ${{ \mathcal E }}_{\mathrm{CR}}$ is the energy deposited in CRs, and ${D}_{{\rm{L}}}$ is the luminosity distance to the source (Waxman & Bahcall 1999; Senno et al. 2018). We consider different values for ${f}_{\mathrm{jets}}$, the fraction of SNe for which the neutrino emission is directed toward the Earth. Depending on the physical model, the emission may e.g., be collimated (relativistic jets; ${f}_{\mathrm{jets}}\ll 1$), or quasi-spherical (shock breakouts; ${f}_{\mathrm{jets}}\sim 1$). We also account for SN misclassifications, in which cases neutrino signals are not injected. For this simplified example, we assume that CC-SNe may be identified from their photometric light curves with a 10% fake-rate (Muthukrishna et al. 2019).

The TS of the stacked signal is computed by averaging the individual neutrino detection significance of all events passing a selection threshold, ${\sigma }_{\mathrm{SN}}^{\min }$. This threshold is defined by the detection significance of the corresponding optical SN. In order to derive the final stacked detection significance, we simulate the background hypothesis; we impose the nominal optical detection procedure, but do not inject any signals into the neutrino data. We generate $\sim {10}^{7}$ such realizations of a full LSST survey, accounting for random coincidence between observables.

The results of the analysis are shown in Figure 3. We compare different values of ${f}_{\mathrm{jets}}$ and ${{ \mathcal E }}_{\mathrm{CR}}$, different optical detection thresholds (${\sigma }_{\mathrm{SN}}^{\min }=3\sigma ,5\sigma$), and different durations of LSST surveys (5, 10 yr). The search would be sensitive in the case of high ${f}_{\mathrm{jets}}$ and ${{ \mathcal E }}_{\mathrm{CR}}$ values, in accordance with Senno et al. (2018). Using our method, it would be possible to e.g., significantly constrain shock breakout scenarios. If no signal is detected, the sensitivity curves could be used to derive physical limits on neutrino emission.

For the neutrino point-source search, the RNN was used to derive joint detection probabilities for two experiments. Similarly in this case, we avoid having to impose a specific relationship between two messengers. In addition, our approach provides a statistical framework to optimize searches, and to derive limits on non-detections. In this illustrative example, we use the RNN to tune the redshift range and to relax the detection threshold of individual events. This circumvents the need to a priori refine the selection criteria of SNe. While we do not explicitly incorporate SN classification or redshift estimation as part of our pipeline, such extensions are feasible (Muthukrishna et al. 2019). They are planned for future work, paving the way for real-time applications.

For SNe engines of sufficient power, relativistic jets manage to break out of the progenitor. Depending on their inclination, these may be observed as long GRBs. Observations of GRBs at high energies are very interesting, e.g., for the study of acceleration mechanisms in relativistic shocks. In many cases, $\gamma$-ray emission extending up to $\mathrm{GeV}$ energies has been detected (Ackermann et al. 2013). Recently, emission of up to hundreds of $\mathrm{TeV}$ has been detected for the first time with ground-based Cherenkov Telescopes, which were following up alerts from other instruments (Abdalla et al. 2019; Acciari et al. 2019).

For the current example, we simulate an uninformed search scenario for the upcoming Cherenkov Telescope Array (CTA; Acharya et al. 2017). We focus on low-luminosity bursts (LL-GRBs), a sub-class of long GRBs (Virgili et al. 2009), which have been connected to SNe having mildly relativistic outflows (Cano et al. 2017), and are potential sources of UHECRs and neutrinos (Murase & Ioka 2013; Boncioli et al. 2019). The true properties of LL-GRBs are not well constrained by observations. Our sample therefore encompasses a wide range of spectral and temporal parameters. As such, it may be used to illustrate the detection prospects of such bursts, rather than to make precise predictions on observable rates. We correspondingly choose simple models, simulating the spectra of LL-GRBs as either simple power laws (PLs) in energy and time, or as PLs having exponential cutoffs.

We simulate observations for the Northern array of CTA using the ctools analysis framework (Knödlseder et al. 2016). The simulations consists of $\gamma$-ray-like events. These correspond to true $\gamma$-rays, as well as to cosmic rays and electrons passing the nominal CTA selection and classification cuts. The inputs to the RNN are event counts, ${n}_{\gamma }$, within $\eta =4$ logarithmically spaced energy bins between 30 and 200 $\mathrm{GeV}$. They are integrated within circular regions of interest (RoIs), over 1 s periods (in order to probe the prompt emission phase of long GRBs). We construct the RNN to represent ${\tau }_{\mathrm{RNN}}=25$ s of data. The first ${\tau }_{\mathrm{enc}}=20$ steps correspond to the background, and the final ${\tau }_{\mathrm{dec}}=5$ steps to the putative signal. In total, the RNN receives $4\times 25=100$ inputs.

For anomaly detection, the RNN is trained to predict $\gamma$-ray-like counts in the absence of signals. In addition, we employ a classification approach. Here, the network is trained with examples of both background and signal events, where signals are injected over the ${\tau }_{\mathrm{dec}}=5$ step interval.

Figure 4(A) shows the real-time discovery potential of the RNN. As expected, bursts with lower redshifts and higher luminosities are more likely to be detected. In general, a large fraction of the parameter space is accessible. In Figure 4(B) we compare our algorithm with the likelihood-based method for source detection of ctools. The RNN performs similarly or better. For anomaly detection, this is achieved without relying on instrument response functions.

Zoom In Zoom Out Reset image size

Contrary to anomaly detection, classification requires some assumption on sources as part of training. However, the performance of the method is shown to be robust to this constraint. In the current example, we simulate the intrinsic spectra of LL-GRBs as exponentially cutoff PLs. The classifier is trained exclusively with simple PL examples. However, nonetheless it is able to generalize, and outperforms the likelihood approach by achieving higher detection rates. We emphasize that for both RNN configurations, training will primarily utilize real CTA data, once it becomes available, rather than simulations. (For instance, for classification a hybrid approach is possible, where simulated signals are injected into background sequences from real data.) Correspondingly, the dependence on instrumental modeling is minimized using our framework.

The advantage of the RNN is particularly evident when only PL models are used as part of the ctools likelihood fit. This represents a realistic strategy for blind searches, for which simple assumptions are generally made. Both our algorithms are comparatively agnostic to the properties of sources. They therefore enable real-time detection of a wider range of transients compared to standard techniques. This illustrates the merits of the methodology for unbiased searches, and the potential for unexpected discoveries.

Our software pipeline and chosen RNN architecture are shown in Figure 1. The RNN is implemented using tensorflow (Abadi et al. 2015). It may be decomposed into two elements, an encoder and a decoder (Cho et al. 2014). The RNN accepts ${\tau }_{\mathrm{RNN}}={\tau }_{\mathrm{enc}}+{\tau }_{\mathrm{dec}}$ time steps as input. The different steps are implemented as LSTM cells. A cell is composed of a pair of LSTM layers, respectively, comprising 128 and 64 hidden units (the set of parameters tuned during training). Each step receives a collection of (analysis-specific) $\eta$ inputs. The inputs are independently normalized, such that their nominal range of values for the background training sample is mapped to the interval, $[0,1]$.

In our nominal approach we employ an anomaly detection technique. We utilize sequences of input data, $S({\tau }_{\mathrm{RNN}},\eta )$, which for training correspond to the response of an instrument in the absence of signal events. The RNN is used to predict the expected background of the experiment, $B({\tau }_{\mathrm{dec}},\eta )$, having the same data structure as $S({\tau }_{\mathrm{RNN}},\eta )$. Transients are then detected as significant divergences from these predictions. Training involves minimizing a loss function, which is defined as the mean squared error between B and S for each set of ${\tau }_{\mathrm{RNN}}$–η.

For classification, an external layer is appended to the output decoder. The latter maps $B({\tau }_{\mathrm{dec}},\eta )$ into logits, $\zeta ({\tau }_{\mathrm{dec}})$. Each of these is a proxy for the RNN probability density function (PDF) for data of a given step to belong to the signal class (Goodfellow et al. 2016). Training proceeds by minimizing a cross-entropy loss function for the different steps, where the final probability is taken as the average, ${\zeta }_{\mathrm{dec}}={\left\langle \zeta \right\rangle }_{{\tau }_{\mathrm{dec}}}$. The corresponding TS is based on the ratio between the background and signal PDFs for a given value of ${\zeta }_{\mathrm{dec}}$ (Cranmer et al. 2015). As for anomaly detection, calibration of test statistics into $p$-values is performed once following the training stage, using simulations.

In the following we detail the data reduction, source modeling, and definition of test statistics used for the example analyses presented in this study. A short summary is given in Figure A1.

Zoom In Zoom Out Reset image size

We search for clustering in two publicly available all-sky $(\mathrm{Decl}.\,\in \,\{-85^\circ ,85^\circ \})$ neutrino samples. The first comprises IceCube event lists of track-like muon neutrino candidates, taken between 2011 and 2012 (MJD 55694–56415), with the IC-86I detector configuration (Aartsen et al. 2017). The second is of ANTARES muon neutrinos, observed within the same time period (Adrián-Martínez et al. 2014). Both samples include time-stamps, reconstructed event directions, and angular uncertainties. In the case of IceCube, energy-proxy metrics and tabulated instrument response functions (IRFs) estimates are available as well. These public data sets are limited in scope. It is therefore not feasible to, e.g., reliably correct for instrumental dead-time, or to estimate the efficiency of event reconstruction and selection cuts. An advantage of our methodology, is that such explicit corrections are not always necessary, so long as unbiased self-consistent observables can be derived.

The nominal inputs to our RNN are based on neutrino event density metrics with respect to a given location on the sky, ${\boldsymbol{x}}$. We define the density for a given neutrino as

$$\begin{eqnarray}&&{s}_{\nu }({\boldsymbol{x}})=\min \left(\displaystyle \frac{{\delta }_{\nu }}{{{\rm{\Delta }}}_{\nu ,{\boldsymbol{x}}}},\,1\right).\end{eqnarray} \tag{ A1 }$$

Here ${\delta }_{\nu }$ is the angular uncertainty of the event, and ${{\rm{\Delta }}}_{\nu ,{\boldsymbol{x}}}$ is the angular distance between the position of the neutrino and ${\boldsymbol{x}}$. We note that this definition of observables, while well motivated, is not unique. However, our framework is designed to provide self-consistent detection probabilities, given such informative input parameters.

We integrate the event densities over 24 hr time periods, which effectively avoids dependence of the event rates on R.A. (Aartsen et al. 2015). For the IceCube sample, the data are split into four logarithmically spaced bins in the energy proxy, ${E}_{\nu }$, between 10 $\mathrm{GeV}$ and 8 $\mathrm{PeV}$. For ANTARES the summation is inclusive in energy. We therefore have five density metrics as RNN inputs,

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{r}\begin{array}{l}{S}_{\nu -{\rm{I}}{\rm{c}}{\rm{e}}}^{\epsilon }({\boldsymbol{x}})=\displaystyle \sum _{1\,{\rm{d}}{\rm{a}}{\rm{y}},\,{\rm{\Delta }}{{\rm{E}}}_{\nu }^{\epsilon }}{s}_{\nu }({\boldsymbol{x}}),\\ \,{\rm{a}}{\rm{n}}{\rm{d}}\,{{\rm{S}}}_{\nu -{\rm{A}}{\rm{N}}{\rm{T}}}({\boldsymbol{x}})=\displaystyle \sum _{1\,{\rm{d}}{\rm{a}}{\rm{y}}}{{\rm{s}}}_{\nu }({\boldsymbol{x}})\end{array},\end{array}\end{array}\end{eqnarray} \tag{ A2 }$$

with energy bins, $\epsilon \,\in \,\{1,2,3,4\}$. The response of the IceCube and ANTARES detectors depends on the zenith of the observed event, which therefore is added as an auxiliary input to the RNN. We nominally use the IceCube events to determine the zenith for a particular sequence. Considering the time of arrival and the location of the observatories, we calculate the corresponding zenith for the same putative source with ANTARES. In cases where an IceCube source is not within the field of view of ANTARES, we use background data instead. In total, we have $\eta =6$ inputs per time step. The collection of inputs is derived for ${\tau }_{\mathrm{RNN}}=15$ days periods. The first ${\tau }_{\mathrm{enc}}=10$ days are assumed to only contain background events. Transient signals are searched for within the next ${\tau }_{\mathrm{dec}}=5$ days interval. In total, the RNN receives $6\times 15=90$ inputs.

We construct a background sample of ${10}^{4}$ ${\tau }_{\mathrm{RNN}}$ sequences, where potential transient signals are removed by scrambling the events in R.A. and time of detection (Aartsen et al. 2015). The RNN is trained to predict the neutrino event densities in each of the five days being probed, for a particular ${\boldsymbol{x}}$. The outputs of the network of the five event densities for a given day are denoted by ${B}_{\nu -\mathrm{Ice}}^{\epsilon (1-4)}$ and ${B}_{\nu -\mathrm{ANT}}$ for the two experiments. We use these to calculate a combined test statistic,

$$\begin{eqnarray}&&{\mathrm{TS}}_{\nu }({\boldsymbol{x}})={\mathrm{TS}}_{\nu -\mathrm{Ice}}({\boldsymbol{x}})+{\mathrm{TS}}_{\nu -\mathrm{ANT}}({\boldsymbol{x}}),\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{TS}}_{\nu -\mathrm{Ice}}({\boldsymbol{x}})=-\displaystyle \sum _{5\,\mathrm{days},\,\epsilon }{\mathrm{log}}_{10}\left(\displaystyle \frac{{B}_{\nu -\mathrm{Ice}}^{\epsilon }}{{S}_{\nu -\mathrm{Ice}}^{\epsilon }}\right),\\ \quad \mathrm{and}\quad {\mathrm{TS}}_{\nu -\mathrm{ANT}}({\boldsymbol{x}})=-\displaystyle \sum _{5\,\mathrm{days}}{\mathrm{log}}_{10}\left(\displaystyle \frac{{{\rm{B}}}_{\nu -\mathrm{ANT}}}{{{\rm{S}}}_{\nu -\mathrm{ANT}}}\right).\end{array}\end{eqnarray} \tag{ A4 }$$

High values of ${\mathrm{TS}}_{\nu }$ correspond to large discrepancies between (background) predictions and the corresponding true data. TS-values are calibrated to $p$-values, ${p}_{\nu }$, as a function of the auxiliary parameter (the zenith, binned into 90 intervals of $2^\circ$ width).

Note that we take a generalizable approach, where the two data sets are combined on equal footing. More complicated schemes are also possible, e.g., by introducing relative weights in the definition of ${\mathrm{TS}}_{\nu }$. Similarly, different timescales may be probed, either by modifying the definition of ${\tau }_{\mathrm{RNN}}$, or that of ${\mathrm{TS}}_{\nu }$.

We also use the outputs of the RNN to perform a correlation test between the IceCube and ANTARES events over the entire period of the data set for each spatial position. For this purpose, we consider ${\mathrm{TS}}_{\nu -\mathrm{Ice}}$ and ${\mathrm{TS}}_{\nu -\mathrm{ANT}}$ individually. As part of the calibration stage, we derive the relation between TS-values and $p$-values for each statistic independently. This is done using simulations (as for ${\mathrm{TS}}_{\nu }$ above).

Our metrics for the correlation analysis are Pearson coefficients, based on the $p$-values of the two samples, ${p}_{\nu -\mathrm{Ice}}$ and ${p}_{\nu -\mathrm{ANT}}$. Explicitly, we define,

$$\begin{eqnarray}&&{\mathrm{TS}}_{\nu -\mathrm{IA}}({\boldsymbol{x}})=\displaystyle \frac{{\sum }_{t}\left({\rho }_{\mathrm{Ice}}^{t}-\left\langle {\rho }_{\mathrm{Ice}}\right\rangle \right)\left({\rho }_{\mathrm{ANT}}^{t}-\left\langle {\rho }_{\mathrm{ANT}}\right\rangle \right)}{\sqrt{{\sum }_{t}{\left({\rho }_{\mathrm{Ice}}^{t}-\left\langle {\rho }_{\mathrm{Ice}}\right\rangle \right)}^{2}}\sqrt{{\sum }_{t}{\left({\rho }_{\mathrm{ANT}}^{t}-\left\langle {\rho }_{\mathrm{ANT}}\right\rangle \right)}^{2}}}.\end{eqnarray} \tag{ A5 }$$

Here ${\rho }_{\mathrm{Ice}}^{t}({\boldsymbol{x}})\equiv -{\mathrm{log}}_{10}{p}_{\nu -\mathrm{Ice}}$, and ${\rho }_{\mathrm{ANT}}^{t}({\boldsymbol{x}})\equiv -{\mathrm{log}}_{10}{p}_{\nu -\mathrm{ANT}}$. The summation is over the entire period of the data set, with $t$ representing a particular 5 days interval for a given sample. We proceed to derive the $p$-value for detection of a correlation signal, ${p}_{\nu -\mathrm{IA}}$. For this purpose, we independently scramble ${\rho }_{\mathrm{Ice}}^{t}$ and ${\rho }_{\mathrm{ANT}}^{t}$ in time and R.A. This procedure is used to create multiple realizations of ${\mathrm{TS}}_{\nu -\mathrm{IA}}({\boldsymbol{x}})$, for which the two samples are uncorrelated. Using these background distributions, we account for spurious correlations, as well as for the number of spatial and temporal grid points.

We simulate a neutrino stacking analysis. The general strategy is to identify an accumulated neutrino over-density (based on ${\mathrm{TS}}_{\nu -\mathrm{Ice}}$) in spatio-temporal coincidence with CC-SNe (Senno et al. 2018). This should be distinguished from the more direct cross-correlation approach, which involves optical follow-up of specific neutrino events (Morgan et al. 2019). A direct analysis is in principle preferred, as it enables the detailed study of a particular source (Aartsen et al. 2018). However, it has two main disadvantages for the case of CC-SNe. First, one must have high confidence that a given neutrino is astrophysical, which generally constrains events to very-high energies. Additionally, the sensitivity is limited by the irreducible contamination of unassociated SNe within the uncertainty region of the neutrino. Given the weak nature of neutrino sources and the large number of potential counterparts, indirect population studies become competitive.

SNe for our stacking analysis are simulated according to projected observations with the upcoming LSST Wide-Fast-Deep survey (LSST Science Collaboration et al. 2009). We utilize the public PLAsTiCC data set, which was created as part of an open data challenge in preparation for LSST (Kessler et al. 2019). The data set represents the projected cadence and observing constraints for LSST. For instance, it includes a prototype scheduler for science program optimization; realistic environmental conditions, such as weather and seeing; maintenance downtime; and instrumental artifacts.

The inputs to the RNN are of two types. For a given time step, the first consists of the IceCube neutrino densities, ${S}_{\nu -\mathrm{Ice}}^{\epsilon }({\boldsymbol{x}})$. The second type consists of LSST signal-to-noise metrics in several optical bands, $b\,\in \,\{g,r,i,z,y\}$. These are defined as

$$\begin{eqnarray}&&{S}_{\mathrm{SN}}^{b}({\boldsymbol{x}})=\displaystyle \frac{{m}_{b}}{\delta {m}_{b}},\end{eqnarray} \tag{ A6 }$$

where ${m}_{b}$ and $\delta {m}_{b}$, respectively, stand for an observed magnitude in band, b, and the corresponding uncertainty. As above, the zenith of observation is added as an auxiliary parameter. We thus have $\eta =10$ inputs per time step in total. We construct ${\tau }_{\mathrm{RNN}}=20$ days data sequences, with the first ${\tau }_{\mathrm{enc}}=10$ days representing the background, and the next ${\tau }_{\mathrm{dec}}=10$ days the search period. In total the RNN receives $10\times 20=200$ inputs.

We begin by creating background samples. These represent periods of joint neutrino and optical observations, in which no transient exists in either messenger. The scrambling procedure described for the neutrino point-source analysis is used to derive ${S}_{\nu -\mathrm{Ice}}$. We extract ${S}_{\mathrm{SN}}$ from simulation periods that lack transient signals, based on the true peak time of SNe light curves. The simulated LSST observations include gaps in observations in some/all bands. We account for missing inputs by randomly interleaving background optical data in such gaps. The background samples are used to train the RNN, using ${10}^{4}$ ${\tau }_{\mathrm{RNN}}$ sequences.

For a particular sky position, we derive individual test statistics for source detection with each messenger. For neutrinos, we use ${\mathrm{TS}}_{\nu -\mathrm{Ice}}({\boldsymbol{x}})$ (here defined over ${\tau }_{\mathrm{dec}}=10$ days intervals). For optical SNe, we have

$$\begin{eqnarray}&&{\mathrm{TS}}_{\mathrm{SN}}({\boldsymbol{x}})=-\sum _{10\,\mathrm{days},\,b}{\mathrm{log}}_{10}\left(\displaystyle \frac{{B}_{\mathrm{SN}}^{b}}{{S}_{\mathrm{SN}}^{b}}\right),\end{eqnarray} \tag{ A7 }$$

where ${B}_{\mathrm{SN}}^{b}$ stands for the output of the network in a given band, b. In addition to high values of ${\mathrm{TS}}_{\mathrm{SN}}$, we impose a constraint on optical detections, that a SN is observed over at least three nights in different bands. The correspondence between ${\mathrm{TS}}_{\nu -\mathrm{Ice}}$ and ${\mathrm{TS}}_{\mathrm{SN}}$ and their respective $p$-values, ${p}_{\nu -\mathrm{Ice}}$ and ${p}_{\mathrm{SN}}$, is derived from simulations as a function of zenith.

In the next step of the analysis, we simulate signal samples. The total number of simulated SNe up to a redshift of 0.06 is $\sim 2000\,{\mathrm{yr}}^{-1}$, following the star formation rate (Bernstein et al. 2012). Of these, about 57% are detected with our RNN with at least $3\sigma$ significance, and 42% with $5\sigma$. We assume that the optically detected SNe are in fact unobserved GRBs. The SNe are generally characterized by the peak time of their light curves, which occurs on average 13 days after the putative GRB (Cano et al. 2017). As the exact time of the emission is uncertain, we randomly generate it as a Poisson process. We consider different values for ${f}_{\mathrm{jets}}$, the fraction of SNe for which the neutrino emission is directed toward the Earth. Depending on the physical model, the emission may e.g., be collimated (relativistic jets; ${f}_{\mathrm{jets}}\ll 1$), or quasi-spherical (shock breakouts; ${f}_{\mathrm{jets}}\sim 1$). We also account for misclassifications of CC-SNe, by introducing a 10% fake-rate, for which neutrino signals are not injected (Muthukrishna et al. 2019).

A neutrino flare is assumed to occur during the short prompt phase of the explosion. We relate the simulated flux of muon neutrinos to the energy deposited by SNe in CRs, ${{ \mathcal E }}_{\mathrm{CR}}$. We use a flat ${\rm{\Lambda }}$ CDM cosmology, with ${{\rm{\Omega }}}_{{\rm{m}}}=0.3$ and the Hubble constant, ${{\rm{H}}}_{0}=70\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. The luminosity distance of optical SNe, ${D}_{{\rm{L}}}$, is derived from their redshift. The fluence of simulated muon neutrinos is then given by ${{ \mathcal F }}_{{\nu }_{\mu }}\propto {{ \mathcal E }}_{\mathrm{CR}}/{D}_{{\rm{L}}}^{2}$ (Waxman & Bahcall 1999; Senno et al. 2018), assuming a PL spectrum for the parent CRs with a spectral index of $-2$, and a flare duration of $\sim 10\mbox{--}100\,{\rm{s}}$. We derive the expected observed neutrino signal from the spectrum using the IRFs of IceCube (Senno et al. 2018).

The test statistic of the stacked signal is defined as

$$\begin{eqnarray}&&{\mathrm{TS}}_{\nu -\mathrm{SN}}({\tau }_{\mathrm{LSST}},{\sigma }_{\mathrm{SN}}^{\min })=\left\langle -{\mathrm{log}}_{10}{p}_{\nu -\mathrm{Ice}}\right\rangle ,\end{eqnarray} \tag{ A8 }$$

where the average is over all SNe having been detected during a survey of duration, ${\tau }_{\mathrm{LSST}}$ (either 5 or 10 yr). In this context, detections are defined as events having an optical detection significance, ${\sigma }_{\mathrm{SN}}$, higher than a given threshold, ${\sigma }_{\mathrm{SN}}^{\min }$ (either $3\sigma$ or $5\sigma$). In order to derive the stacked detection $p$-value, ${p}_{\nu -\mathrm{SN}}$, we simulate the background hypothesis; we impose the nominal optical detection procedure, but do not inject any signals into the neutrino data. We generate $\sim {10}^{7}$ such realizations of a full LSST survey of duration ${\tau }_{\mathrm{LSST}}$, accounting e.g., for misidentified SNe, and for random coincidence between observables (Senno et al. 2018).

Observations are simulated for the Northern array of CTA using ctools (Knödlseder et al. 2016), one of the proposed analysis frameworks for the observatory. The simulations produce $\gamma$-ray-like events. These correspond to true $\gamma$-rays, as well as to CRs and electrons, which pass all selection and classification cuts. ctools allows generation of background (CR, electron) sky-maps, as well as of background+signal observations, by inclusion of source spectral models. We use IRFs optimized for 30 minutes observations at zenith angles of $20^\circ$. The RoI for the simulation is chosen as a circular area with a radius of 0.25°. It is centered at the position of the putative source, and is displaced by 0.5° from the center of the camera.

The inputs to the RNN are $\gamma$-ray-like event counts, ${n}_{\gamma }$, within $\eta =4$ logarithmically spaced energy bins, E_γ, between 30 and 200 $\mathrm{GeV}$ ($\epsilon \,\in \,\{1,2,3,4\}$), integrated within 1 s intervals within the RoI,

$$\begin{eqnarray}&&{S}_{\gamma }^{\epsilon }=\displaystyle \sum _{1\,{\rm{s}},\,{\rm{\Delta }}{E}_{\gamma }^{\epsilon }}{n}_{\gamma }.\end{eqnarray} \tag{ A9 }$$

For anomaly detection, we train the network with background-only ctools simulations. The RNN therefore predicts ${B}_{\gamma }^{\epsilon (1-4)}$, the background event counts in each of the energy bins for each of the steps. We construct the RNN to represent ${\tau }_{\mathrm{RNN}}=25$ s of data. The first ${\tau }_{\mathrm{enc}}=20$ steps correspond to the background, and the final ${\tau }_{\mathrm{dec}}=5$ steps to the putative signal. In total, the RNN receives $4\times 25=100$ inputs. We train the network using ${10}^{4}$ background sequences.

Our framework allows for explicit statistical models to be used, avoiding the need for calibration with background simulations. We illustrate this for the LL-GRB analysis with Poissonian models of the signal and background+signal hypotheses. These are defined as ${P}_{\mathrm{pois}}(k| \lambda )={e}^{-\lambda }{\lambda }^{k}/k!$, given an observed number of events, k, and a rate, $\lambda$. We approximate $\lambda$ from the event count corresponding to a given hypothesis, integrating over the final ${\tau }_{\mathrm{dec}}=5$ steps of the RNN for a given energy bin,

$$\begin{eqnarray}&&{\lambda }_{S}^{\epsilon }=\sum _{5\,{\rm{s}}}{S}_{\gamma }^{\epsilon },\quad \mathrm{and}\quad {\lambda }_{B}^{\epsilon }=\sum _{5\,{\rm{s}}}{{\rm{B}}}_{\gamma }^{\epsilon }.\end{eqnarray} \tag{ A10 }$$

We then define the test statistic,

$$\begin{eqnarray}&&{\mathrm{TS}}_{\gamma -{\rm{A}}}=-2\sum _{\epsilon }\mathrm{ln}\left(\displaystyle \frac{{P}_{\mathrm{pois}}({S}_{\gamma }^{\epsilon }| {\lambda }_{B}^{\epsilon })}{{P}_{\mathrm{pois}}({S}_{\gamma }^{\epsilon }| {\lambda }_{S}^{\epsilon })}\right).\end{eqnarray} \tag{ A11 }$$

We directly derive the pre-trials $p$-values for anomaly detection, ${p}_{\gamma -{\rm{A}}}$, using Wilks's theorem (Wilks 1938), which was validated by background simulations).

In addition to anomaly detection, we employ a classification approach. In this case, we require examples of signal events for training. To simulate the possible $\gamma$-ray signatures of LL-GRBs for CTA, we nominally model their spectra as simple PLs in energy, E, and time, t,

$$\begin{eqnarray}&&{M}_{\mathrm{PL}}(E,t)\propto {E}^{-{\rm{\Gamma }}}{t}^{-T}.\end{eqnarray} \tag{ A12 }$$

The spectral and temporal indices, Γ and T, are parameters of the models. The true properties of these events are not well constrained, due to the scarcity of observations. We therefore consider a wide range of parameters. Motivated by bright bursts, detected at high energies with Fermi-LAT, we randomize over likely values, $1.9\lt {\rm{\Gamma }}\lt 2.7$ and $0.8\,\lt T\,\lt 2$ (Ackermann et al. 2013). The flux normalization is randomly shifted with regard to these reference events in redshift and luminosity to the expected ranges for LL-GRBs (Acharya et al. 2017). The observed spectra are corrected for interactions with the extragalactic background light (EBL), which attenuates high-energy $\gamma$-rays (Franceschini et al. 2008).

The network is trained for classification using ${10}^{4}$ background sequences and ${10}^{4}$ signal sequences. Signals are injected over the ${\tau }_{\mathrm{dec}}=5$ step interval, as part of the ctools simulation. The network is trained with a wide range of such models having different flux normalizations, corresponding to different signal-to-noise ratios with respect to the background. The signal models incorporate PL temporal decay. Correspondingly, late-time models are equivalent to early-time models with relatively lower flux normalization. The inclusive composition of the training sample therefore enhances generalization (time-invariance) of the RNN.

While the inputs to the RNN for classification are the same as for anomaly detection, the output in this case is a single number, ${\zeta }_{\mathrm{dec}}$. The latter takes low values for input sequences that do not include transient signals, and high values in the presence of high-significance signals. We use the training sample to derive ${P}_{\mathrm{CLAS}}^{{\rm{B}}}$ and ${P}_{\mathrm{CLAS}}^{{\rm{S}}}$, the approximated PDFs for the background and for the background+signal hypotheses (see Figure A2 (A)). Under the assumption that ${\zeta }_{\mathrm{dec}}$ is monotonic with the ratio of PDFs (Cranmer et al. 2015), the classification test statistic can be defined as

$$\begin{eqnarray}&&{\mathrm{TS}}_{\gamma -{\rm{C}}}({\zeta }_{\mathrm{dec}})=-2\mathrm{ln}\left(\displaystyle \frac{{P}_{\mathrm{CLAS}}^{{\rm{B}}}({\zeta }_{\mathrm{dec}})}{{P}_{\mathrm{CLAS}}^{{\rm{B}}}({\zeta }_{\mathrm{dec}})+{P}_{\mathrm{CLAS}}^{{\rm{S}}}({\zeta }_{\mathrm{dec}})}\right).\end{eqnarray} \tag{ A13 }$$

We calibrate the relation between ${\mathrm{TS}}_{\gamma -{\rm{C}}}$ and the corresponding pre-trials ${p}_{\gamma -{\rm{C}}}$ with background and signal simulations, using Wilks's theorem. The calibration is defined as a function of ${\zeta }_{\mathrm{dec}}$, as illustrated in Figure A2(B).

Zoom In Zoom Out Reset image size

In a realistic scenario, an uninformed search would be performed for RoIs at different positions across the field of view of CTA. The search would be repeated multiple times, depending on the definition of ${\tau }_{\mathrm{dec}}$ and on the amount of observing time. Correspondingly, the detection significance must be corrected for trials. We take, ${p}_{\mathrm{post}}=1\,-{\left(1-{p}_{\mathrm{pre}}\right)}^{{n}_{\mathrm{trial}}}$, as the relation between pre- and post-trials probabilities. We use the number of trials, ${n}_{\mathrm{trial}}\sim 3.6\times {10}^{7}$, corresponding to searches over 100 hr of observations, over the entire CTA field of view.

It is currently difficult to estimate the rate of false detections for CTA, as the observatory will explore a new regime of sensitivity (Acharya et al. 2017). Experience with real data will improve the fidelity of detections. As a first step, we performed systematic checks on the stability of the $\gamma$-ray-like event counts we use as input to the RNN. The counts are susceptible to fluctuation due to imperfect $\gamma$-ray reconstruction, and to uncertainties on the IRFs. In particular, energy dispersion below $\sim 50\,\mathrm{GeV}$ may result in migration between bins, and can change the energy threshold of the analysis. We studied these effects by comparing simulations where we vary the IRFs by their expected uncertainty (of up to $10 \%$). We found that the propagated uncertainties on counts do not significantly affect our results. We also tested theoretical uncertainties on the effect of the EBL. Comparing different models of the EBL (Franceschini et al. 2008; Dominguez et al. 2011; Gilmore et al. 2012), we found negligible impact on the observed spectra of LL-GRBs. This is primarily due to the low redshift and energy regimes we consider in the current study.

We also simulate a separate signal sample, where LL-GRBs are modeled as PLs having exponential cutoffs,

$$\begin{eqnarray}&&{M}_{\mathrm{EC}}(E,t)={M}_{\mathrm{PL}}\cdot \exp \left(-\displaystyle \frac{E}{{E}_{\mathrm{cut}}}\right).\end{eqnarray} \tag{ A14 }$$

In this case, we scan a range of cutoff energies, $1\lt {E}_{\mathrm{cut}}\,\lt 120\,\mathrm{GeV}$. The cutoff models are not used for training. Rather, we utilize them to illustrate the robustness of our methods (see Figure 4). As discussed above, the networks perform well for different source types. For instance, they enable identification of GRBs having cutoff spectral models, despite only being trained with simple PL examples.