Combined detection of supernova neutrino signals

Supernova neutrino detection is usually implemented as a real-time trigger system to provide an immediate reaction to the occurrence of signal-like structures in a stream of data. This real-time detection regime implies continuous computation of the probability that the observed data at some point in time is consistent with the expected neutrino signal coming from the late stages of the evolution of massive stars: either the burst of neutrinos from a core collapse, or the cooling of a nearby red supergiant core in the Si-burning presupernova phase.

The most thorough approach for locating the signal in the presence of background would be fitting the measured data with the supernova signal parameters (distance, mass, starting time etc.). However, this approach has limited efficacity because existing simulations might not cover all the physics possibilities and hypotheses. Moreover, performing fits to multiple models in a large parameter space is computationally intensive and appears not feasible for the case of the real-time supernova detection, as such evaluation has to be done repeatedly for each new portion of data received. Therefore we consider a more simple and robust procedure of hypotheses test, measuring the supernova significance as the deviation from the background and avoiding the scan of all supernova signal parameters, except for the starting time $t^{*}$.

Consider an experiment that measures timestamps of individual neutrino interaction candidates $\{t_{i}\}$. The task of the supernova detection system is to tell if the observed data contains a signal from a supernova and to evaluate the starting time of that astrophysical event. Assuming the supernova starting time $t^{\ast}$, such system scans the incoming data and performs statistical tests to estimate the significance of the signal starting at $t^{\ast}$ in the considered data slice. To discriminate between two hypotheses, $H_{0}$ (background only) and $H_{1}$ (background + supernova), one has to define a test statistic as a function of the incoming data: $\ell(\{t_{i}\})$.

The significance of observing deviations from the null hypothesis $H_{0}$ using the test statistic $\ell$ is defined by the tail of the distribution:

$$p(\ell)=\int\limits_{\ell}^{\infty}P(\ell^{\prime}|H_{0})\,d\ell^{\prime}.$$

(2.1)

For convenience this can be converted to a $z$-score, defined as

$$z(\ell)=\Phi^{-1}(1-p(\ell)),$$

(2.2)

where $\Phi$ is the cumulative standard normal distribution function. This z-score is equivalent to that derived in Gaussian statistics as the number of sigma away from the mean.

The supernova triggering system would react on the signal present in the data when the significance z-score exceeds a defined threshold $z_{thr}$. Thus the efficiency of the supernova detection under hypothesis $H$ is defined by the probability distribution above this threshold:

$$\varepsilon(H)=\int\limits_{z_{thr}}^{\infty}P(z|H)\;dz.$$

(2.3)

Calculation of the significance values (2.1) and (2.2) requires knowing the probability distribution of the test statistic under a given hypothesis $P(\ell|H)$.

For a hypothesis $H$ predicting the event rate $r(t)$, the number of events $n$ within the considered time window $[t^{*},t^{*}+\Delta t]$ is a Poisson random number around the expected rate:

$$P(n|H)=\frac{e^{-R}R^{n}}{n!},\;\text{ where }R=\int\limits_{t^{*}}^{t^{*}+\Delta t}r(t)\,dt.$$

(2.6)

Test statistics in both CA and SA are additive functions of the timestamps of neutrino candidate data $\{t_{i}\}$: $\ell(\{t_{i}\})=\sum_{i}\ell(t_{i})$. Thus the total test statistic distribution can be expressed as a weighted sum of contributions for each $n$:

$$P(\ell|H)=\sum_{n=0}^{\infty}P(\ell|n,H)\cdot P(n|H).$$

(2.7)

Assuming that all the individual events are causally independent, we can express $P(\ell|n,H)$ as a convolution with a single event distribution $P(\ell|1,H)$:

$$P(\ell|n,H)=P(\ell|n-1,H)*P(\ell|1,H),$$

(2.8)

which, using the convolution theorem, becomes:

$$P(\ell|n,H)=\mathcal{F}^{-1}\left\{\Phi^{n}(z)\right\},$$

(2.9)

where $\Phi(z)=\mathcal{F}\left\{P(\ell|1,H)\right\}$ is a characteristic function for $P(\ell|1,H)$, and $\mathcal{F}\left\{\right\}$ and $\mathcal{F}^{-1}\left\{\right\}$ denote direct and inverse Fourier transformations.

The distribution of the test statistic value for a single event with timestamp $t$ is defined as:

$$P(\ell|1,H)=\int\limits_{t^{*}}^{t^{*}+\Delta t}dtP(t|1,H)\int\limits_{0}^{\infty}d\ell^{\prime}\delta(\ell^{\prime}-\ell(t)).$$

(2.10)

and the distribution of a single timestamp within the time window is defined by the event rate, predicted by hypothesis $H$: $P(t|1,H)=r(t)/R$. These expressions are used for the numerical computation of the characteristic function $\Phi(z)$.

Putting (2.6) and (2.9) into (2.7) we get:

$$P(\ell|H)=\sum\limits_{n=0}^{\infty}\frac{e^{-R}R^{n}}{n!}\mathcal{F}^{-1}\left\{\Phi^{n}(z)\right\}=\mathcal{F}^{-1}\left\{e^{-R}\sum\limits_{n=0}^{\infty}\frac{(R\Phi(z))^{n}}{n!}\right\}.$$

(2.11)

Finally, combining the series into to an exponent, we get the expression for the test statistic distribution:

$$P(\ell|H)=\mathcal{F}^{-1}\left\{e^{R[\Phi(z)-1]}\right\}.$$

(2.12)

In case of multiple experiments each using their own test statistic functions $\{\ell_{n}(t)\}$, their combination is nontrivial. However, if each of $\ell$ is a log likelihood ratio $\ell(t)=\log\frac{P(t|H_{1})}{P(t|H_{0})}$ and all the experiments use the same hypotheses $H_{0}$ and $H_{1}$, these values are additive:

$$\ell_{comb}=\sum\limits_{n=1}^{N_{exp}}\ell_{n}(\{t_{i}^{n}\})=\sum\limits_{n=1}^{N_{exp}}\sum_{i}\ell_{n}(t_{i}^{n}),$$

where $N_{exp}$ is number of experiments to combine.

In this case the distribution of combined test statistic $P(\ell_{comb}|H)$ becomes a convolution of distributions for the individual experiments, so (2.12) becomes:

$$P(\ell_{comb}|H)=\mathcal{F}^{-1}\left\{\prod_{n=1}^{N_{exp}}\exp\left(R_{n}[\Phi_{n}(z)-1]\right)\right\}.$$

(2.13)

The shape analysis method described above is implemented as a python package and is available at [7].

Briefly, the computation of the resulting significance requires the following preparation steps:

1. 1.

   Define the background $B$ and expected signal $S$ rates.

2. 2.

   Compute the single event distribution $P(\ell|1,H_{0})$ using (2.10), and its image $\Phi(z)$ with the fast Fourier transformation.

3. 3.

   Calculate the test statistic distribution for the background-only hypothesis $P(\ell|H_{0})$ using eqs. (2.12) or (2.13) in case of many experiments.

Then the data analysis can be done, based on the events timestamps $\{t_{i}\}$ for each considered supernova starting time $t^{*}$:

1. 4.

   Calculate the test statistic value $\ell(t^{*},\{t_{i}\})$ via (2.5)

2. 5.

   Evaluate the significance corresponding to the obtained value $\ell$ using $P(\ell|H_{0})$ via (2.1), (2.2)

This results in a time series of significance values, representing the likelihood that a supernova following the $H_{1}$ hypothesis started at each point in time.

Computationally the shape analysis method is more complicated than the counting analysis, which relies only on the Poisson distribution. However, the calculation of the test statistic distribution needs to be done only when the background level (or expected signal) has changed, so in practice these reevaluations can be done with large time intervals, leading to a smaller computational load and latency on average. For example, in the case of the NOvA supernova triggering system [5] the reevaluation of the background level and recalculation of the null hypothesis distribution is performed every 10 minutes, as that is the characteristic time over which background rates in NOvA are stable.

Figure 1 shows an example of such computation for an experiment with a background rate $B=0.1\,\mathrm{event/s}$ and a signal with total expected $\int S(t)dt=3\,\mathrm{events}$, injected at time $T=0$. For the signal shape we used an analytical expression to imitate a typical core-collapse anti-neutrino signal:

$$S(t)\sim(1-e^{-t/\tau_{0}})\cdot e^{-t/\tau_{1}},$$

(2.14)

with time parameters $\tau_{0}=0.1\,\mathrm{s}$, $\tau_{1}=1\,\mathrm{s}$.

The length of the time window for both analyses is set to 5 s, so each neutrino interaction (vertical lines in the top panel) increases the test statistic value for the assumed supernova start time within 5 s in the past. For the counting analysis this impact is constant within the whole time window, while for the shape analysis it is defined by the signal and background shapes, thus giving a more precise prediction of the start time, and a higher maximal significance.

A comparison of the shape analysis performance with the standard counting analysis is made by calculating the significance for picking out the supernova signal above noise at various distances (figure 3, right). The band shows the 68 % uncertainty assuming Poisson fluctuations of the background and signal event numbers.

The profile of the supernova significance vs. the assumed time of the SN start (figure 3, left) shows that SA has maximum significance when the assumed time matches the actual start of the signal, while CA has a less pronounced peak, spread to as far as the counting analysis window (in case of low statistics, like in NearDet). This, and the symmetrical shape of the significance profile in the SA case allows a more precise determination of the SN burst start.

These advantages of SA over CA remain true (with a slight decrease in total significance) even in the case when expected signal model is different from the detected signal shape (as shown in [5]). This is caused by the fact that different models still share the general features of the core-collapse signal.

The performance of the shape analysis method depends on the choice of the hypotheses $H_{0}$ and $H_{1}$, which are defined by the background $B(t)$ and signal $S(t)$ rates in (2.5). While background can be studied from the detector data (as in NOvA case), the signal shape is obtained from the supernova simulations, introducing a model uncertainty.

Here we consider four signal models introduced in table 1 as possible shapes of the actual supernova signal $S_{true}(t)$.

We use (2.12) and (2.1 – 2.3) to calculate the distribution of the significance value $z$ for the hypothesis, defined by an event rate $r(t)=B(t)+S_{true}(t)$, while using the same background $B(t)$ and selecting an expected signal model $S(t)$ for the definition of the test statistic function (2.5). This allows us to probe the analysis performance for various combinations of expected and input signal models $S(t)$ and $S_{true}(t)$.

Table 2 summarizes the mean expected significance for the supernova signals at $d=10\,\mathrm{kpc}$ and the distance corresponding to $50\%$ efficiency, using the shape analysis with various signal hypotheses $S(t)$ (the four supernova models, “simple” parametrization), and counting analysis, vs. simulated input signal models $S_{true}(t)$.

The table shows that the best results in terms of the significance and detection distance are achieved when the received signal shape matches the expected one. However, even the “simple” parametrization reflecting only the most basic features of the core-collapse supernova signal shape provides a good improvement in comparison to the counting analysis.

One of the tasks of the supernova detection system is the reconstruction of the supernova signal start time $t^{*}_{rec}$ for the given set of data $\{t_{i}\}$.

As the significance $z(t^{*},\{t_{i}\})$ measures a goodness of fitting the data with the supernova signal starting at the time $t^{*}$, we use the time of maximal significance as an estimator for $t^{*}$:

$$t_{rec}^{*}=\operatornamewithlimits{argmax}_{t^{*}}z(t^{*},\{t_{i}\}).$$

(3.1)

In order to study the precision of such estimator, we perform a large number of toy Monte-Carlo simulations of the neutrino event timestamps $\{t_{i}\}$, with each neutrino time drawn from a distribution following the event rate $R(t)=B+S_{true}(t)$, where $S_{true}(t)$ is one of the input signals from $d=10\,\mathrm{kpc}$ distant supernova starting at a given time $t^{*}_{true}$. Then we analyze these simulated datasets assuming various signal hypotheses in the same way as in previous section, calculating the $t^{*}_{rec}-t^{*}_{true}$. Figure 4 shows such distributions for the simulated signal from $9.6\,\mathrm{M_{\odot}}$ NH signal in the NOvA far detector, processed by the three different analyses. Since the tails of the distribution of the $t_{rec}$ defined in (3.1) depend on the considered time window, one cannot use the standard deviation or the distribution quantiles in order to evaluate the distribution width (and thus $t_{rec}^{*}$ estimator precision). Instead we use the time of the distribution maximum and the full width on the half maximum (FWHM) of the distribution peak as a metric of the method precision.

The time precision estimation procedure has been applied to all combinations of expected and received signals as in the layout of table 2. Typical FWHM values of the obtained distributions accord with the example depicted in figure 4 ($\sim 1000\,\mathrm{ms}$ for CA and $\sim 200\,\mathrm{ms}$ for SA) showing a robust $\sim 5$ times improvement in timing precision when accounting for the signal shape.

Currently the triggering systems on NOvA detectors work independently of each other. However combining signals from NearDet and FarDet can improve the total sensitivity of the system, as shown in figure 5.

Possible time shifts between the detectors might affect the combined significance, however most experiments use GPS time synchronization with the uncertainty below 100 ns[10], which is much lower than the scale of the significance peaks on figure 5. Moreover, since the detection is based on pattern matching rather than determination of the sharp rise of the signal front, these shifts also have negligible impact to the precision of the determination of $t^{*}_{rec}$.

Table 3 summarizes the results of the calculation of the expected NOvA reach, assuming detection threshold $z_{thr}=5\,\mathrm{\sigma}$ and detection efficiency $\varepsilon=50\,\mathrm{\%}$, for various signal models and depending on the applied analysis. This calculation shows that using the shape analysis for individual detectors allows detection about 1 kpc further than counting analysis. An additional increase of about 1 kpc is achieved by applying the combined shape analysis.

As in the case of the core-collapse supernova neutrino burst, the presupernova neutrinos can be detected by various neutrino and dark matter experiments in several interaction channels [14]. Here we will limit our consideration to inverse beta-decay (IBD) channel in three experiments: KamLAND, Borexino and SK-Gd.

The detection efficiencies of the considered detectors vs. anti-neutrino energy $E_{\nu}$, as well as presupernova anti-neutrino energy spectra in the considered models, are depicted in figure 6. Figure 7 shows the resulting presupernova event rates versus time expected in the KamLAND experiment for the models considered in this work.

Table 4 summarizes the signal event rates expected for different presupernova models, background levels, and the time windows chosen for the counting analysis, for each experiment considered.

To demonstrate features of the shape analysis applied to the presupernova neutrino search, we hereafter consider the Odrzywolek model and normal neutrino hierarchy, unless noted otherwise.

While the detection of the core collapse supernova neutrino signal would provide the precise timing and other physical information of the supernova explosion process, the presupernova signal has an additional goal of predicting the supernova event beforehand. Thus the metric of how far in advance we can detect the supernova is as important as the maximal detection distance.

The shape analysis probes how well the signal+background shape fits the acquired data, so the optimal discrimination (maximal significance) is reached when the full-length signal is present in the data. However, since the signal shapes in figure 7 show almost monotonic behaviour in time, even an incomplete early part of the signal in the data can provide a considerable significance elevation when fitted with the full signal shape.

In other words, applying the shape analysis improves significance even at the early stages of the presupernova signal. Figure 8 demonstrates that the shape analysis allows an earlier detection of the presupernova signal from 200 pc distance than the counting analysis. The panels in figure 9 show the dependency of the maximal presupernova significance (1 minute prior to the collapse) on the supernova distance, for the same model.

Figure 10 shows the regions where the detection efficiency is above 50 %, 90 % and 99 % in the counting and shape analyses for the Odrzywolek model with NH, assuming the detection threshold $z_{thr}=5\,\mathrm{\sigma}$. The exact value of this threshold in real detection system should be defined by the tolerable false alarm probability.

Joint analysis of the detectors using  (2.13) provides an additional improvement of the presupernova signal detection efficiency. Figure 11 shows the regions with 90 % efficiency for separate shape analyses and their combinations. Including of the analyses with low sensitivity to the presupernova signal, such as Borexino or SK-Gd neutron mode, still increases the combined reach. For example, the 90 % efficiency of detection of a presupernova signal at 200 pc distance occurs about two hours earlier for Borexino+KamLAND, than only for KamLAND. However the combination doesn’t provide large benefits at the very early stages ($>6\,\mathrm{hours}$ before a supernova), because the shape analysis is optimized on fitting the full signal shape.

The accuracy and performance of the shape analysis method depends on the knowledge of the background and signal rates. Since the actual signal is not known in advance, we study the uncertainty of the method by calculating the sensitivity reaches for various signal models and compare them to the counting analysis. Figure 12 shows that while the maximal reach is always obtained when the expected signal and the received signal coincide, other models still show better reach than the counting analysis. This reflects the fact that even a wrong signal model describes the presupernova signal better than a constant rate in a given time window (as in the counting analysis).

Table 5 presents quantitative metrics of the expected reach, prediction time and significance using various analysis methods for the detection modes in considered detectors and their combinations. The performance was calculated for all above-mentioned presupernova models and neutrino mass hierarchies using the Odrzywolek NH model as the expected signal in SA. Distance and time in this table correspond to the 50 % efficiency detection reach (in contrast to 90 % in figures 11, 12).