Identification of neutrino bursts associated to supernovae with Real-time Test Statistic (RTS$^2$) method

This paper proposes a new approach for the selection of low-energy neutrino bursts, such as the ones detected after a supernova. It exploits the temporal structure of the expected signal with respect to the more diffuse background by defining a"Real-time Test Statistic"(RTS) that would allow identifying very weak signals, hard to select using standard clustering methods. For a given background rate, the new method (RTS$^2$: RTS for Supernovae) increases signal efficiency while keeping the same false alarm rate for Poisson-distributed background. By adding a spatial penalty term to the definition of RTS, one can also reject spatially-correlated backgrounds such as the ones due to spallation events. Furthermore, the method is easy to implement in a real-time monitoring system as RTS can be computed recursively for successive events, and it can be easily adapted for detectors of all scales that may want to send prompt alerts e.g. through SNEWS 2.0 network.


I. INTRODUCTION
The detection of low-energy neutrino bursts is particularly relevant in the search for supernovae. Prompt detection would allow providing early warning for the subsequent electromagnetic observations, that is crucial to understand the underlying mechanisms.
The usual method consists in watching for upper fluctuations from the typical expected background. In most of the neutrino detectors, the main background at energies relevant for supernovae can be modelled with Poisson statistics, i.e. a constant rate with simple statistical fluctuations. For a given background rate r, the number of background events in a time window w is simply following the Poisson distribution with λ = r × w: Poisson(m; r, w) = e −r×w (r × w) m m! .
Let's denote τ phys the typical time scale of the physical phenomenon under study. For instance, τ phys ∼ 10 s is expected in the case of supernovae. If rτ phys 1, Gaussian approximation can be used and the search for neutrino bursts can be simplified to a search of an excess with a z-score [2].
However, if the detector has a relatively low background (rτ phys 1), this method cannot be applied and the usual technique consists in dividing the data in clusters using time windows of size w, covering the physical time scale τ phys . For supernovae, w = 20 s is commonly used. There are two different approaches to define these time windows, as illustrated in Figure 1: • Sliding windows of width w, with each window starting in the middle of the previous one [3].
• Dynamic windows of width w, one starting from each selected event [1,8] It is then straightforward to define a cut on the cluster multiplicity in order to achieve a target false alarm rate FAR (e.g. 1/year or 1/century) by using the Poisson survival function.
For sliding time windows, and for dynamic time windows, where the factor in front of the sum is the "rate" of time windows. For instance, for a background rate of 0.01 s −1 and a false alarm rate of 1 century −1 , one gets multiplicity cuts of 7 and 8 events for the sliding and dynamic methods respectively (16 and 17 for r = 0.1 s −1 ). Nevertheless, these methods are limited as neither take into account that signal bursts will have a very different time structure (localised in a few seconds interval) than the background (expected to be uniformly distributed within the time window).
Recently, there have been attempts in using additional time-related variables, such as the time difference ∆t between the first and the last event in a cluster [6]: typically ∆t ∼ τ phys for signal-only cluster and ∆t ∼ w arXiv:2103.09733v1 [astro-ph.HE] 17 Mar 2021 for background-only cluster. This can however suffer from background contamination, in particular at high background rates.
In this publication, an alternative approach, that is getting rid of the traditional time window definition while benefiting from the time structure of signal burst, is proposed. The section II will present the new method. In section III, the results when applied to toy SN simulation will be presented. In section IV, it is shown that the method can still be applied in case of remaining non-Poisson background, with slight modifications. The section V presents a brief discussion about the different aspects of the new method.

II. METHOD
Each selected neutrino event i is characterised by a contribution f (t; t i ), maximal at t = t i . One can then compute: where the sum is performed over all events. This continuous function F can then be used to search for an excess i.e. F(t) > F cut . The value of the cut can be optimised to obtain a given false alarm rate. If N centuries of background are simulated and the requested false alarm rate is one per century, one should find F cut such that F(t) exceeds this threshold less than N times in total.
If one uses this is fully equivalent to the dynamic time window approach presented in section I. However, in this publication, the following function will be considered: where T c is a characteristic time. F(t) is re-written as: For the implementation in an online process, one may simply store in memory the value F(t i ) at the previous event and rewrite F for t i < t < t i+1 : and: Average value: In the case of Poisson-distributed background, the average value of F(t) is: where n is the number of events within the period [0, t] and P (t i ) is the distribution of the time of event i. Using order statistic and assuming the events are uniformly distributed between 0 and t, the latter may be written as: By taking the limit t → ∞ (n → rt), one can derive the full computation and obtain: This ensures that F remains stable in case of background-only events.
For simplicity, in the following, the method is named RTS 2 for Real-time Test Statistic for Supernovae.
The Figure 2 shows the typical distribution of F(t) (using T c = 10 s) if it is computed regularly for background-only scenario (with r = 0.05 s −1 ), compared with the distribution of the maximum value reached for a supernova signal (see section III A for the description of the latter). As expected from the Equation 12, the background distribution keeps relatively small values and the average value is F ∼ 0.50 = rT c .

III. RESULTS
In this section, various background rates will be considered. The performance of the new method presented in section II will be compared to ones of the simple clusterisation in dynamic time windows (Equation 3), referred as the "multiplicity" method.
In the following, T c = 10 s will be used as it is found to be a good compromise in the search for supernova signal with similar time scale.

A. Signal simulation
The simulation is performed in a generic way, so that it is relatively independent on the supernova models and on the neutrino detector: • A given number of signal events n S is supposed to be detected. This could be converted to a supernova distance for a given detector (fixed mass, target and efficiency) using SNOwGLoBES [5] ; • Their time distribution is parametrised by a double-exponential shape as in [6]: with typically τ 1 = 10-100 ms and τ 2 1 s, representing respectively the rising time and the decaying time of the signal. In the following, when not explicitly specified, τ 1 = 10 ms and τ 2 = 1 s are used, which fit well the current SN models and SN1987a data.
Each signal burst is simulated on a short time interval (typically τ = 30-50 s), on top of simulated background events with rate r bkg . The signal efficiency for a given cut C is obtained by simulating a large number of signal bursts and by computing the fraction of bursts satisfying C during the period τ (one cluster passing the multiplicity cut or F(t) exceeding given threshold).

B. Raw performance
Independently on the signal simulation, one can first compare directly the cuts applied on F and on multiplicity in order to achieve a given false alarm rate respectively in the RTS 2 and "multiplicity" methods, as illustrated in Figure 3. This is relevant as these cuts correspond to the minimal number of events needed to trigger an alert.
The figure clearly shows that the new method is reaching lower multiplicities and is more robust with increasing multiplicities. The penultimate point on the right in Figure 3 corresponds to the scenario presented in the Figure 2, with r = 0.05 s −1 and with a typical cut F 4 − 5 as already visible in the histogram. Cut to achieve FAR

C. Signal selection efficiency
The Figure 4 presents the obtained signal efficiency for different fixed background rates and number of injected signal neutrinos. It confirms the results presented in the previous paragraph: the RTS 2 method performs much better for low-multiplicity signal bursts and the effect is visible for all background rates. Referring back to the example of r = 0.05 s −1 from the Figure 2, one can achieve non-negligible efficiencies for n sig > 6, as seen in the clear separation between background and signal in the histogram. For a fixed false alarm rate and background rate, one may then reduce the detection threshold (and therefore increase the distance horizon) of the supernova search.

D. Stability with signal shape
The Figure 5 illustrates how the signal efficiency changes with different signal shapes (Equation 13) and different sizes of time window for the "multiplicity" method. Reducing the time window down to 5 s allows increasing the efficiency with respect to the longer duration, but it does not reach the performance of the newly proposed method.

IV. NON-POISSON BACKGROUND
Neutrino detectors may suffer from spatially correlated backgrounds whose times will not be distributed with a Poisson distribution. This is the case of spallation events, where radioactive decays may follow the interaction of high-energy cosmic ray muons with the matter of the detector. In Super-Kamiokande, this is handled by checking the spatial distribution of the identified vertices after clustering [1], computing the "dimensionality" of the vertex distribution and only selecting clusters identified as "volume-distributed" (vertices uniformly distributed in the volume of the detector).
One may try to adapt the real-time test statistic defined in section II in order to penalise events if they are detected nearby previous events, so that spatially-correlated clusters are effectively suppressed: where p i is the weight applied to event i. A naive approach is to define: where { x i } are event positions, T SP and σ SP are the time and spatial characteristic of the spatial penalty and the sum runs over the events with t j < t i (it may be sufficient to only consider events with t i − t j 10T SP ). The formula ensures that p i ∼ 1 if there are no events nearby event i, while p i → 0 otherwise (event will not contribute to F). The weight of a given event is reduced even further if more than one other event is detected nearby.
For a correlated cluster of N events and characterised by a multi-normal distribution with σ = σ B , one can derive that F SP (t c ) N i=1 p i ≤ 2, if σ B = σ SP and assuming all events are arriving at the same time t c , even for N → ∞. The effect of the cluster is effectively suppressed. The term σ SP can easily be tuned to the scale of the background of interest in a given detector.
The Figure 6 illustrates how the implementation of the spatial penalty terms allows to reduce the contamination due to non-Poisson background, as peaks related to non-Poisson background (point-like sources) are shrunk, while signal bursts remain at the same level than without spatial penalty. The Figure 7 presents the rejection power of RTS 2 method for spatially-correlated clusters.
One can see that, even though the chosen spatial penalty term is relatively simple, it is quite effective to reject point-source background. It is only in the case of low false alarm rates (e.g., 1/hour) that extra care may be needed. Furthermore, the definition of p i = p({x j } j≤i , {t j } j≤i ) may be optimised to reject specific backgrounds related to a given experiment, as the proposed formula only takes care of point-like sources while e.g. spallation background may have vertices distributed along a line.
To classify each newly detected neutrino as part of a background "burst" (p i → 0) or not (p i ∼ 1), the use of a k-nearest neighbours algorithm or other outliers detection techniques may be very promising.

V. DISCUSSION
The presented approach would allow to select supernovae bursts with relatively low multiplicity (∼ 5 − 10), effectively increasing the horizon compared to more standard approaches.   For reference, based on numbers from Table 4 of [4], the Hyper-Kamiokande detector would be able to reach M31-Andromeda (d ∼ 778 kpc) while smaller current detectors would still trigger for most of the Milky Way satellite galaxies [7] (d 200 kpc).
The approach can be easily be implemented in a realtime monitoring system, as a recursive process of the type F(t i+1 ) = a(t i+1 − t i )F(t i ) + p i+1 , where a(τ ) = e −τ /Tc characterises the decay of the signal and the {p i } characterise the weights of individual events (either = 1 or < 1 if spatial penalties are applied).
The threshold may be adjusted to safely reduce contamination by using Monte Carlo simulations of the background or by data-driven methods. The latter are particularly relevant if one wants to achieve false alarm rates of 1/(day, week, month), as limited datasets would be sufficient to validate it.

VI. CONCLUSION
In this paper, the new RTS 2 method has been presented and applied to the detection of neutrino bursts with a generic time spectrum that fit well current supernovae models. The performances were compared to using solely a cut on the multiplicity in fixed size time window and it has been shown that the new method performs well, especially for higher background rates.
The method allows high flexibility on the choice of the function characterising the contribution of each event and it may be further improved by considering adding new variables in the analysis, such as the event energy.
It is possible to simply adapt the RTS 2 formula in order to take care of the non-Poisson distributed backgrounds, either with a basic penalty term (e.g. Equation 15) or with specially adapted methods, to reject specific backgrounds in a given detector (outlier detection techniques or machine learning approaches).
Not only could it help to extend the horizon of existing and future large neutrino detectors beyond their current reach, but it may also be useful for smaller detectors that need a consistent way to define and send alerts to the community, in particular through SNEWS 2.0 [4]. The latter aims to gather data from neutrino detectors all over the world (with various sizes and targets) for combination and initiate multi-messenger follow-up of a candidate supernova. This would ensure taking the most from the next galactic core-collapse supernovae that the astronomy community is waiting for. ACKNOWLEDGMENTS I would like to thank the local group in the Department of Physics of the University of Padova, as well as the astrophysics working group of the Hyper-Kamiokande collaboration, for useful discussions. CloudVeneto is acknowledged for the use of computing and storage facilities. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 754496.